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ABSTRACT 


Described  is  the  behavior  of  a  rotorcraft  equipped  with  Higher  Harmonic  Stall 
Control  (HHSC)  and  a  Reverse  Velocity  Rotor  (RVR).  Current  rotorcraft  are  limited  in 
forward  flight  speed  by  retreating  blade  stall  and  compressibility  effects  on  the  advancing 
blade.  Stall  occurs  as  the  blade  encounters  increasingly  severe  reverse  flow.  HHSC 
enables  conventional  rotor  systems  to  fly  on  the  forward  and  aft  sections  of  the  rotor  disk, 
greatly  reducing  reliance  on  the  mixed  flow  regions  defined  by  the  advancing  and 
retreating  blades.  Employment  of  the  RVR  allows  lift  generation  while  the  rotor  is 
experiencing  reverse  flow.  A  similar  type  of  two  per  revolution  (2/rev)  input  can  be 
tailored  to  deliver  maximum  benefit  to  RVR  equipped  rotorcraft. 

Modification  of  the  Joint  Army  Navy  Rotorcraft  Analysis  and  Design  (JANRAD) 
computer  program  allows  2/rev  cyclic  input,  use  of  the  RVR,  and  analysis  using  high 
fidelity  graphical  output  to  examine  angle  of  attack,  coefficient  of  lift,  and  air  load. 
Computational  results  show  performance  gains  in  conventional  helicopters  and  high 
speed  flight  potential  for  RVR  equipped  aircraft.  The  RVR  is  applied  to  the  Joint  Heavy 
Vertical  Lift  (JVHL)  aircraft  conceptual  design  for  preliminary  analysis.  This  conceptual 
design  can  be  used  as  an  indicator  of  the  performance  of  a  high  speed  RVR  equipped 
aircraft. 
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I.  INTRODUCTION 


A,  MOTIVATION  FOR  RESEARCH 

The  next  generation  of  ship-based  military  vertieal  heavy  lift  rotoreraft  must  be 
eapable  of  eondueting  high  speed  long-range  eombat  assaults.  It  must  be  able  to  operate 
in  all  weather  conditions,  day  or  night,  from  prepared  or  unprepared  surfaces. 
Operational  requirements  are  shown  in  the  figures  below. 

Long  Range  Heavy  Lift  Aircraft  Mission  Profile 


Figure  1 .  Mission  Profile. 


Vehicle  Gross  Weights  for  Air  Transport 
70  I - , - ^ ^ ^ - , - 


HMMWVLW155  M198  M7VR  LAV  HEMAT  AAAV  M1A1 
vehicle 


Figure  2.  Payload  Requirements. 
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These  requirements  are  outside  the  envelope  of  any  existing  Vertieal/Short  Take- 
Off  and  Landing  (VSTOL)  aireraft  to  date.  By  the  year  2020,  the  CH-46E  will  have  been 
replaced  by  the  MV-22A.  But,  the  MV-22A  does  not  have  sufficient  internal  or  external 
payload  capacity  to  transport  heavy  equipment  or  large  quantities  of  bulk  logistics 
support.  The  CH-53E  does  not  have  sufficient  internal  load  capacity  to  carry  heavy 
equipment  such  as  the  Mk28/Mk928  7-ton  truck  needed  to  support  USMC  field  artillery 
ashore.  It  does  not  have  sufficient  range  with  large  external  loads  like  the  Eight  Armored 
Vehicle  (EAV)  to  support  the  Operational  Maneuver  from  the  Sea  (OMETS)  and  Ship  to 
Objective  Maneuver  (STOM)  operations  up  to  200  NM  inland  envisioned  in  the  future. 
Although  a  modified  CH-53X  has  been  proposed,  no  manufacturing  plan  has  been 
formalized  and  current  plans  call  for  the  phase  out  of  existing  CH-53Es  by  2025. 
Clearly,  a  new  aircraft  must  be  developed  to  meet  future  requirements. 

B.  HIGHER  HARMONIC  STALL  CONTROL  AND  REVERSE  VELOCITY 

ROTORS 

Two  of  the  most  compelling  technologies  that  might  be  used  to  meet  these 
operational  requirements  are  Higher  Harmonic  Stall  Control  (HHSC)  and  Reverse 
Velocity  Rotors  (RVR).  Both  of  these  technologies  have  been  studied  in  great  detail  but 
have  never  enjoyed  widespread  acceptance  within  the  helicopter  industry.  An  aircraft 
using  HHSC  and  the  RVR  system  would  retain  all  the  characteristics  of  a  conventional 
helicopter  at  low  airspeeds  but  would  also  be  capable  of  high  speed  flight  (in  excess  of 
200  knots). 

HHSC  differs  from  Higher  Harmonic  Control  (HHC)  in  the  following  manner. 
HHC  is  an  active  control  concept  whose  principal  objective  has  been  to  reduce  helicopter 
airframe  vibrations.  More  recently  it  has  been  applied  to  improve  helicopter  performance 
as  well.  Eor  HHC  principle  frequencies  of  interest  are  (n  -1)Q  ,  nCl ,  and  {n  +  \)Q.  in  the 
rotating  system  and  nO.  in  the  fixed  fuselage  system.  Here,  n  ,  is  the  number  of  blades 
and  fl  is  the  rotor  rotational  speed.  [Wood  et.  al,  1985]  In  contrast  to  HHC,  HHSC 
relates  usually  to  2Q  or  second  harmonic  frequency.  [Arcidiacono,  1961]  More  recently 
it  has  been  extended  to  3Q  frequency  as  well.  The  objective  of  HHSC  when  applied  to 
the  rotor  in  high-speed  flight  is  to  smooth  out  the  distribution  of  lift  and  eliminate  stall 
over  the  entire  rotor  disk. 
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A  preliminary  design  of  an  RVR  based  vertical  heavy  lift  aircraft  was  developed 
but  analysis  and  verification  of  rotor  performance  was  nearly  impossible  since  traditional 
design  methodology  and  performance  calculations  worked  well  with  the  rotorcraft  only 
up  to  170  knots.  The  planned  cruise  speed  of  the  RVR  based  aircraft  was  300  knots. 
After  several  attempts  to  constrain  our  design  to  work  with  traditional  rotor  design  tools  it 
was  determined  the  best  course  of  action  would  be  to  develop  a  new  way  to  investigate 
the  behavior  of  the  rotor  system  in  which  both  RVR  and  HHSC  technology  would  be 
applied. 

C.  JANRAD  2P  RVR 

Due  to  the  unique  nature  of  RVR  and  HHSC  the  original  Joint  Army  Navy 
Rotorcraft  Analysis  and  Design  (JANRAD)  program  was  modified  to  allow  for  HHSC, 
use  of  RVRs,  and  several  other  enhancements.  The  core  pitch-thrust  moment  iteration 
technology  that  forms  the  core  of  JANRAD  was  unchanged. 


Figure  3.  Graphical  Depiction  of  Pitch- Thrust  Moment  Iteration.  [From: 

Gerstenberger  and  Wood,  1963] 

A  complete  description  of  the  original  program,  assumptions,  and  operating 
constraints  can  be  found  in  theses  by  Nicholson  [Nicholson,  1993]  and  Eccles.  [Eccles, 
1995]  Sections  two  and  three  discuss  the  employment  and  effects  of  using  the  HHSC  and 
RVR  systems.  These  chapters  will  provide  a  working  knowledge  of  how  each  system 
works  by  outlining  program  modifications,  assumptions,  and  actual  JANRAD  2P_RVR 
output.  UH-60  aircraft  data  will  be  used  to  further  illustrate  applicability  to  current 
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helicopter  design.  JANRAD  2P_RVR  was  used  for  the  analysis  of  the  Joint  Vertical 
Heavy  Lift  Aircraft  (JVHL)  outlined  in  section  four.  Complete  user  instructions  for 
JANRAD  2P_RVR  are  provided  in  section  five. 

The  study  concludes  with  a  brief  list  of  recommendations  for  improvement  and 
suggestions  for  further  study.  A  complete  listing  of  all  the  modules  required  to  run  the 
JANRAD  2P_RVR  program  is  given  in  the  appendices. 
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II.  HIGHER  HARMONIC  STALL  CONTROL 


A.  INTRODUCTION 

Since  the  early  1960s  researchers  have  pursued  the  idea  of  ineorporating  two  per 
revolution  (2/rev  or  2P)  eontrol  inputs  to  rotor  systems.  Studies  ranged  from  inereasing 
forward  fight  speed  to  the  reduetion  of  aireraft  vibrations.  Computer  eontrolled  Higher 
Harmonie  Control  (HHC),  3P,  4P  and  5P,  showed  signifieant  vibration  reduetion  and  was 
the  subjeet  of  a  detailed  flight  test  program  eondueted  by  Hughes  Helieopters,  Ine.  under 
a  NASA/ Army  eontract.  [Wood  et.  al,  1985]  For  the  study  of  this  thesis  the  emphasis 
will  be  on  performance  airspeed  inereases  and  eontrolling  the  rotor  stall  pattern  in  whieh 
2P  and  3P  Higher  Harmonie  Stall  Control  (HHSC)  will  be  applied. 

B,  BLADE  PITCH  MOTION  AND  THE  SWASHPLATE 

Blade  piteh  motion  eomes  from  two  sourees,  namely  eommanded  eolleetive  and 
eyelie  piteh  input  from  the  helieopter  eontrol  system  and  elastie  deformations  (twist)  of 
the  blade  and  eontrol  system.  This  section  will  focus  on  the  cyclie  eommanded  input. 
Elastie  deformations  are  negleeted  at  this  level  of  study. 

The  piteh  time  history  of  the  rotor  blade  ean  be  represented  by  a  Fourier  series  as 
follows: 

6  {y/)  =  6*0  +  6^^  oos(^^)  +  6^^  sin(^)  +  ...  +  6^^  cos{ni//)  +  0^^  + ... 

where  0  is  the  blade’s  feathering  angle,  y/  is  the  blade  azimuth  position,  and  y/  =  Q.t. 
Conventional  eontrol  inputs  produeed  by  the  swashplate  eonsist  of  eolleetive  piteh,  6*o , 
and  the  first  harmonies  of  the  Fourier  series:  the  lateral  eyelie  0^^  and  the  longitudinal 
eyelie  0^^ .  A  conventional  blade  piteh  (feathering)  motion  ean  be  deseribed  by  the 
Fourier  series 

0  (y/)  =  00+  01^  eos(^)  +  0^^  sin(^)  =  0o  +  0^^  cos(Qt)  +  sin(Qt) 

This  one  per  revolution  eontrol  (1/rev  or  IP)  formed  the  eore  of  JANRAD’s 
original  trim  routine.  The  0^^  term  eontrols  the  lateral  orientation  of  the  rotor  disk  and 

6*;^  term  eontrols  the  longitudinal  orientation.  In  the  above  equation,  t  =  0  has  been 
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taken  with  the  blade  in  the  trailing  or  aft  ( ^  =  0  )  position.  Remember  that  for  a  rotor 
with  a  centrally  located  flap  hinge  there  is  an  exact  90  degree  force/displacement  phase 
lag.  In  the  case  of  pure  6^^  (sine)  longitudinal  cyclic  motion,  the  minimum  aerodynamic 
flapping  motion  occurs  at  ^  =  90  degrees,  where  the  minimum  flapping  displacement 
occurs  90  degrees  later  at  ^  =  180  degrees;  the  maximum  aerodynamic  flapping  moment 
occurs  at  y/  =  270  degrees,  and  the  maximum  flapping  displacement  occurs  90  degrees 
later,  at  =  360  degrees. 

The  swashplate  is  the  key  to  affecting  the  pitch  control  of  the  rotor  blades.  The 
swashplate  has  rotating  and  nonrotating  disks  concentric  with  the  rotor  shaft.  Both  disks 
slide  up  and  down  in  response  to  collective  inputs  or  they  can  be  tilted  to  an  arbitrary 
orientation  in  response  to  cyclic  inputs  made  by  the  pilot.  The  pilot  cyclically  changes 
the  pitch  of  the  blades  about  the  swashplate.  HHSC  inputs  can  be  incorporated  into 
existing  helicopters  by  modification  of  the  swashplate.  HHSC  inputs  would  be 
transmitted  to  the  rotor  using  a  stationary  outer  member  with  a  track  attached  to  a 
conventional  swashplate  assembly.  This  differs  dramatically  from  HHC  systems  which 
rely  on  computer  controlled  actuators  connected  to  the  swashplate  to  alter  the  rotor  pitch. 
Figure  4  compares  a  HHC  system  with  a  proposed  HHSC  mechanism. 


Figure  4.  HHC  Actuator  System  [From:  Wood  et.  al,  1985]  Compared  to  a  HHSC 
Swashplate  and  Ring  System.  [From:  Arcidiacono,  1961] 
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C.  EXPERIMENTAL  ANALYSIS  USING  JANRAD  2P  RVR 

Since  research  by  Arcidiacono  has  suggested  a  significant  increase  in  maximum 
forward  speed  of  a  helicopter  could  be  obtained  through  the  use  of  HHSC  control  of  rotor 
pitch  the  first  modifications  to  JANRAD  were  made  to  the  Fourier  series  that  defines 
cyclic  pitch.  First  iterations  of  JANRAD  2P_RVR  used  the  following  equation  for 
0  given  by 

e  {y/)  =  6»o  +  cos{y/)  +  6,^  siniy/)  +  6^^  cos{2y/)  +  0^^  ^m{2y/)  +  0^^  sin(3^/^  + 

This  equation  incorporates  lateral  and  longitudinal  2P  signals  and  longitudinal  3P 
signals  plus  a  phase,  0^^^^^ .  The  2P  and  3P  0  coefficients  are  user  defined  fixed  values 

while  the  IP  coefficients  are  determined  by  JANRAD  based  on  first  harmonic  trimming 
rules.  Since  the  objective  of  our  study  was  to  maximize  performance  it  was  found  that 
the  best  equation  for  0  was  found  to  be 

0  (y/)  =  00+  cos(^^)  +  sin(^^)  +  0^^  cosi2y/)  +  0,^  sm{3y/  + 

The  2P  sine  term  was  eliminated  because  it  was  causing  increased  rotor  stall  on 
the  left  and  right  side  of  the  rotor  disk,  in  the  exact  location  were  we  wished  to  reduce 
rotor  stall.  The  3P  sine  term  produced  striking  results  and  was  retained  for  experimental 
value,  but  in  the  majority  of  cases  this  term  was  assigned  a  value  of  zero. 

D,  TEST  USING  PROUTY’S  HELICOPTER 

As  stated  earlier,  the  primary  reason  for  using  HHSC  (2P)  is  to  improve  the  lifting 
efficiency  of  the  rotor  by  modifying  the  stall  pattern.  Before  attempting  to  analyze  2P 
behavior  with  JANRAD  2P_RVR,  conventional  first  harmonic  inputs  were  verified  by 
comparing  basic  output  with  known  results.  Prouty’s  sample  helicopter  was  used  as  the 
baseline  for  comparison.  To  determine  if  JANRAD  2P_RVR  was  producing  correct 
output,  rotor  disk  angle  of  attack  contour  plots  were  compared  to  known  solutions.  These 
plots  were  chosen  because  the  stall  patterns  are  directly  linked  to  the  angle  of  attack 
distribution.  All  2P  and  3P  terms  were  set  to  zero  for  the  test.  Figure  5  shows  the 
calculated  first  harmonic  cyclic  pitch  of  the  Prouty  example  helicopter. 
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Figure  5.  Rotor  Azimuth  versus  Harmonic  Cyclic  Pitch  for  the  Prouty  Example 

Helicopter. 

Figure  6  compares  the  angle  of  attack  contour  plot  from  Prouty’s  text  to  JANRAD 
1P2P3P  output. 


«  -  iw 


Figure  6.  Prouty’s  Example  Helicopter  Contour  Plot  (Feft)  Compared  to  JANRAD 
2P_RVR  Output  (Right)  Using  Prouty’s  Example  Helicopter. 

The  similarities  between  the  two  plots  indicate  JANRAD  2P_RVR  output  is 
stable  and  in  close  agreement  with  Prouty’s  numerical  results.  Next,  2P  and  3P  signals 
were  introduced.  Figures  7  and  8  show  the  result  for  a  2P  signal  and  3P  signal, 
respectively. 
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Figure  7.  2P  Signal  Effects  on  Cyclic  Pitch  (Left)  and  Angle  of  Attack  (Right). 


Figure  8.  3P  Signal  Effects  on  Cyclic  Pitch  (Left)  and  Angle  of  Attack  (Right). 


Three  dimensional  lift  distribution  plots  were  used  to  further  verily  the  effects  of 
these  higher  harmonic  signals.  The  plots  are  compared  to  conventional  lift  distributions 
in  Figures  9  and  10. 


Figure  9.  Conventional  (IP)  Lift  Distribution  (Left)  Compared  to  2P  Lift 

Distribution. 
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These  plots  clearly  show  how  the  2P  signal  moves  the  aerodynamic  forces  to  the 
forward  and  aft  sections  of  the  rotor  disk.  By  doing  this  the  helicopter  can  ‘fly’  on  the 
forward  and  aft  sections  of  the  disk  and  therefore  not  rely  as  heavily  on  the  left  and  right 
sections.  The  rotor  is  no  longer  trying  to  provide  propulsive  and  lifting  force  from  the 
unstable  left  and  right  rotor  regions.  This  translates  into  the  ability  to  gain  forward 
velocity.  As  was  mentioned  earlier,  the  3P  signal  was  ineffective  in  showing  any 
performance  gains  but  was  retained  for  experimental  value. 

E.  COMPARISON  TO  ARCIDIACONO 

One  additional  test  was  run  to  verify  JANRAD  2P_RVR’s  integrity.  Thrust 
distribution  and  angle  and  angle  of  attack  distribution  were  compared  to  ideal  thrust  and 
angle  of  attack  distributions  derived  by  Arcidiacono.  [Arcidiacono,  1961]  Figures  1 1  and 
13  show  Arcidiacono’s  original  plots  and  his  attempts  to  meet  the  optimum  distribution. 
Note  that  Arcidiacono’s  sign  conventional  is  opposite  to  the  sign  conventional  used 
throughout  this  study. 
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THRUST  COEFFICIENT,  ONE  BLADE 


NO  HIGHER  HARMONIC  FEATHERING 


AZIMUTH  ANGLE  ,  .  OEG 


Figure  1 1 .  Arcidiacono’s  Results  with  No  2P. 


115.1492  Kts.  advance  ratiQ=0. 29772.  2P  input=0  degs.  TPP  angle=0. 090769  degs.  3P  input=0  degs.  3P  offset=0  degs 


Rotor  Tip  Path  Plane  Stall  Pattern,  advance  ratio=0. 29772 


X  Actual  and  Ideal  Thrust  Distribution  at  Stall  advance  ratio=0. 29772 


Figure  12.  Results  Using  JANRAD  2P_RVR  with  No  2P. 
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^2^  =  -IDEG,  Bgg  =-4  DEG,  Aj^  =  +  IDEG,  Bjg  =-0.5  DEG 


AZIMUTH  ANGLE, ,  DEG 


A2g=-IDEG,  B2g=-4DEG,  A3g=+IDEG,  B3g=-0.5DEG 


0  9  0  180  270  360 

AZIMUTH  ANGLE,  i^.OEG 


Figure  13.  Arcidiacono’s  Results  with  2P  and  3P.  [From;  Arcidiacono,  1961] 


115-1492  Kts.  advance  ratio=0-29776-  2P  inpLit=4  degs.  TPP  angle=0-089514  degs.  3P  input=0  degs.  3P  offeet=0  degs 
Rotor  Tip  Path  Plane  Stall  Pattern,  advance  ratio=0.29776 


X  10"^  Actual  and  Ideal  Thrust  Distribution  at  Stall,  advance  ratio=0.29776 


Figure  14.  Results  Using  JANRAD  2P_RVR  with  2P. 


Figure  14  shows  that  JANRAD  results  were  favorable  in  terms  of  thrust 
coefficient.  Angle  of  attack  distribution  was  expected  to  be  less  favorable  due  to  the 
constraints  imposed  by  JANRAD  for  tip  losses.  Overall,  JANRAD  2P_RVR  validated 
itself  for  use  with  2P  inputs. 
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F.  FORWARD  FLIGHT  VELOCITY  INCREASES 

Previous  versions  of  JANRAD  had  been  validated  for  use  with  the  UH-60A; 
therefore  it  was  selected  as  the  helicopter  to  test  for  potential  speed  increases.  [Eccles, 
1995]  All  UH-60A  parameters  were  loaded  into  JANRAD  2P_RVR.  Sikorsky’s 
proprietary  airfoil,  the  SC1095R8,  was  used  as  the  rotor  blade.  Testing  was  confined  to 
conventional  helicopter  configurations  only. 

In  short,  an  increase  in  maximum  flight  speed  of  14  knots  was  achieved  using  a 
four  degree  fixed  2P  input.  This  is  just  short  of  a  ten  percent  improvement.  No 
significant  difference  in  maximum  forward  speed  was  observed  with  2P  inputs  greater 
than  4  degrees.  This  increase  was  observed  under  the  1 .5  percent  iterative  deviation  rule 
imposed  by  JANRAD  2P_RVR.  That  is,  collective  and  cyclic  pitch  settings  are  adjusted 
until  the  thrust  vector  is  within  1.5  percent  of  the  previous  iteration.  This  stringent 
requirement  ensures  realistic  simulation  within  the  trim  module. 

The  maximum  airspeed  without  any  2P  input  was  159  knots.  The  rotor  is  under  a 
great  deal  of  stress  to  produce  all  the  required  forces  for  flight.  Figures  15  through  18 
depict  the  rotor  behavior  under  these  conditions.  Any  increase  in  speed  caused  JANRAD 
2P_RVR  to  diverge  and  return  to  the  user  a  ‘WILL  NOT  TRIM’  message. 
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Figure  16.  UFI-60A  Lift  Distribution  at  159  Knots. 
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Figure  17.  UFI-60A  Angle  of  Attack  Distribution  at  159  Knots. 
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Figure  18.  UH-60A  Cl  Distribution  at  159  Knots. 

The  maximum  airspeed  with  the  four  degree  2P  input  was  173  knots.  The 
familiar  2P  pattern  is  elearly  evident  in  Figure  21.  The  other  figures  show  results  similar 
in  appearanee  to  the  plots  derived  from  Prouty’s  example  helieopter.  Any  inerease  in 
speed  eaused  JANRAD  2P_RVR  to  diverge  and  return  to  the  user  a  ‘WILL  NOT  TRIM’ 
message. 
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Figure  19.  UH-60A  2P  Cyclic  Pitch  at  173  Knots. 


Lift  Distribution  At  173.2245  Kts,  advance  ratio=0. 38913,  2P  input=4  degs,  3P  input=0  degs,  3P  offsBt=0  degs 
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Figure  20.  UFI-60A  2P  Lift  Distribution  at  173  Knots. 
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Figure  21 .  UFI-60A  2P  Angle  of  Attaek  Distribution  at  173  Knots. 


Figure  22.  UFI-60A  2P  Cl  Distribution  at  173  Knots. 

Although  not  shown,  testing  was  also  carried  out  with  compound  helicopters.  Use 
of  2P  systems  with  compound  helicopters  equipped  with  auxiliary  thrust  devices 
performed  very  similarly  to  conventional  2P  equipped  helicopters  while  the  rotor  is 
providing  forward  thrusting  power.  As  the  auxiliary  thrust  system  begins  to  provide  the 
majority  of  propulsive  force  at  airspeeds  above  180  knots,  if  no  power  was  supplied  to 
the  rotor,  the  disk  would  assume  an  aft  tilting  tip  path  plane  and  rotor  drag  would 
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increase  significantly.  This  situation  is  avoided  completely  in  actual  compound 
helicopters  in  high-speed  flight  by  providing  a  relatively  small  amount  of  power  to 
stabilize  the  rotor,  while  keeping  the  tip  path  plane  at  a  near  zero  angle  of  attack. 

The  application  of  second  harmonic  control  is  not  a  new  idea  for  improving  the 
speed  performance  of  helicopters.  This  thesis  has  verified  with  JANRAD  2P_RVR  that  a 
2P  rotor  system  does  trim  at  airspeeds  up  to  10%  higher  than  conventional  IP  rotor 
systems.  The  next  section  explores  the  use  of  Reverse  Velocity  Rotors  coupled  with  a 
HHSC  control  signal. 
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III.  EMPLOYMENT  OF  REVERSE  VELOCITY  ROTORS 


Current  helicopters  are  limited  in  their  speed  by  three  effects,  namely  stalling  of 
the  tip  of  the  retreating  blade,  loss  of  lift  due  to  reverse  flow  over  the  inboard  portion  of 
the  retreating  blade,  and  compressibility  on  the  tip  of  the  advancing  blade.  [Ewans, 
1973]  The  severity  of  these  effects  is  a  function  of  the  rotational  speed  of  the  rotor 
blades  and  the  forward  velocity  of  the  aircraft.  Note  that  if  we  treat  retreating  blade  stall 
by  slowing  the  rotor  we  also  reduce  compressibility  effect  on  the  advancing  blade. 
Several  studies  have  been  conducted  to  determine  the  best  way  to  combat  retreating  blade 
stall.  Everything  from  exotic  blade  shapes,  segmented  rotors  [Zientek,  2001]  and  forced 
air/circulation  control  have  been  attempted. 

A,  REVERSE  VELOCITY  ROTORS 

A  relatively  new  and  promising  concept  recently  introduced  is  called  the  Reverse 
Velocity  Rotor  (RVR).  RVR  systems  are  a  recent  development,  invented  by  Harold 
Eemont  to  address  the  aerodynamic  issues  contributing  to  retreating  blade  stall  and 
advancing  blade  compressibility  phenomena,  the  primary  factors  that  currently  limit 
helicopter  maximum  forward  flight  speeds.  Instead  of  eliminating  the  reverse  flow  that 
occurs  on  the  retreating  side  of  the  blade  the  RVR  is  actually  designed  to  operate  in  the 
reverse  flow  region.  With  conventional  rotor  blades,  this  would  completely  eliminate  the 
potential  to  generate  lift.  In  fact,  a  ‘reversible’  airfoil  section  with  a  rounded  trailing  edge 
is  employed  to  give  reasonable  lift  and  drag  characteristics  in  the  reverse  flow.  Reverse 
velocity  rotor  systems  are  built  around  double-ended  airfoils  as  shown  in  Eigure  23. 


Typical  RVR  airfoil 


Eigure  23.  RVR  Cross  Section.  [Erom;  Ashby,  2002] 

Reverse  velocity  rotor  systems  represent  a  revolutionary  approach  to  high-speed 
Vertical  Take  Off  and  Eanding  (VTOE)  configurations.  Application  of  this  concept 
makes  possible  rotorcraft  designs  that  are  capable  of  attaining  speeds  significantly  greater 
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than  conventional  or  even  compound  helicopters.  Initial  analysis  shows  cruise  speeds  in 
excess  of  300  knots  are  possible  with  RVR  systems.  Sikorsky  Aircraft  Corporation  has 
completed  a  study  that  shows  future  Runway  Independent  Aircraft  (RIA)  equipped  with 
RVR  systems  could  attain  cruise  speeds  over  300  knots,  double  the  cruise  speed  of 
conventional  helicopters.  [Ashby,  2002]  This  airfoil  design  minimizes  the  impact  of 
retreating  blade  apparent  velocity  reduction  caused  by  summing  the  blades  rotational 
velocity  and  the  aircraft’s  forward  flight  velocity  as  illustrated  in  Figure  24  below. 


Vapparent  =  275  +  507  =  782  ft/S 


Figure  24.  Rotor  Disk  Apparent  Velocities.  [From;  Ashby,  2002] 

As  might  be  expected  the  RVR  system  does  require  somewhat  more  power  to 
hover  and  fly  at  low  airspeeds.  For  this  reason  the  best  application  for  RVR  systems  is 
for  rotorcraft  that  will  spend  the  majority  of  time  in  cruise  configuration,  i.e.  civilian  or 
military  transport.  Aircraft  that  will  spend  the  majority  of  their  flight  time  in  low  speed 
flight  or  hovering,  i.e.  rescue  aircraft,  should  continue  to  employ  conventional  rotor 
blades.  Figure  25  shows  predicted  power  requirements. 
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Power  (hp) 


Figure  25.  Power  Requirements  versus  Velocity.  [From;  Zientek,  2001] 


B,  WIND  TUNNEL  TEST 

Fairchild  Republic  Division  conducted  wind  tunnel  tests  at  NASA  Ames  in  1973 
of  a  1/7  scale  12%  thickness  double  ended  airfoil  at  equivalent  airspeeds  of  over  310  kts. 
[Ewans,  1973]  The  test  rig  is  shown  in  Figure  26. 


Figure  26.  Installation  of  RVR  Test  Rig  in  the  NASA  Ames  12  Foot  Wind  Tunnel. 

[From;  Ashby,  2002] 
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Testing  was  considered  a  success  although  some  large  unbalanced  forces  were 
noted  at  the  edge  of  the  test  envelope.  Data  collection  focused  on  measuring  rotor 
performance,  particularly  maximum  lift  and  effective  lift/drag  ratio  over  the  full 
envelope.  [Ewans,  1975]  The  most  important  data  for  this  study  are  given  in  Figures  27 
and  28.  From  these  plots  data  tables  were  generated  which  correlated  Mach  number, 
angle  of  attack,  and  Cl  or  Cd.  These  tables  were  inserted  into  a  MATFAB  function 
which  used  a  two-dimensional  interpolation  to  provide  the  JANRAD  2P_RVR  thrust 
moment  and  drag  moment  routines  with  the  appropriate  Cl  and  Cd  values. 


SECTION  urr  coErnciENT  throcoh  ss<r  angle  or  attack 

is%  AZRrOIL  -  FULL  SCALE  REYNOLDS  NUMBER 


UR  CoMndwit 


Figure  27.  RVR  Cl  versus  Angle  of  Attack.  [From;  Ewans,  1975] 
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SECTION  DRAG  COEFFICIENT  TOROUGH  360*  ANGLE  OF  ATTACK 
12%  AIRFOIL  -  MODEL  REYNOLDS  NUMBER 


Figure  28.  RVR  Cd  versus  Angle  of  Attack.  [From;  Ewans,  1975] 

C.  EXPERIMENTAL  ANALYSIS  USING  JANRAD  2P  RVR 

JANRAD,  in  its  original  form,  was  not  set  up  to  calculate  trim  parameters  for  a 
rotorcraft  traveling  at  300  knots.  Several  new  assumptions  had  to  be  made  to  make 
JANRAD  2P_RVR  compatible  with  high  speed  RVR  systems. 

1,  Advance  Ratio 

Before  outlining  the  assumptions,  the  advance  ratio  parameter  must  be  defined. 
The  advance  ratio  is  the  ratio  of  the  rotor’s  horizontal  speed  to  its  tip  speed.  The  tip 
speed  ratio,  /u ,  is  defined  as; 

H  =  cos  a  /  QR 

where 

=  forward  airspeed 

a  =  tip  path  plane  angle 
Q  =  rotational  velocity  in  radians  per  second 
R  =  rotor  radius 

With  a  conventional  rotor  system  as  jj,  exceeds  0.4,  collective  and  cyclic  pitch 
values  were  unable  to  satisfy  rotor  lift,  rotor  drag,  or  resultant  thrust  vector  requirements. 
The  system  falls  out  of  ‘trim’  with  respect  to  pitching,  rolling  and  yawing  moments.  Even 
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after  applying  second  or  HHSC  the  program  was  unable  to  trim  above  advance  ratios  of 
.5.  This  is  not  surprising  when  we  recognize  that  at  //  =  0.5,  one  half  of  the  retreating 
blade  {y/  =  270  degrees)  is  in  reverse  flow.  In  contrast,  RVR  systems  operate  in  an 
environment  with  advance  ratios  ranging  between  1.0  and  2.0.  This  fivefold  increase  in 
advance  ratio  in  RVR  aircraft  over  aircraft  equipped  with  conventional  rotors  illustrates 
the  significant  differences  between  the  two  systems. 

2,  RVR  Assumptions  and  Requirements 

Several  assumptions  were  made  to  enable  computational  trimming  of  an  RVR 
system  using  the  traditional  JANRAD  trimming  program.  For  high-speed  rotorcraft 
flight  (//  >  0.5)  JANRAD  2P_RVR  only  allows  RVR  trimming  on  compound  helicopters 
where  a  wing  provides  the  majority  of  vertical  lift  at  cruise  speeds.  Proper  employment 
of  a  wing  can  reduce  rotor  loads  up  to  80  percent.  The  reduction  in  rotor  loading 
eliminates  blade  stalling  concerns  due  to  high  loading.  It  is  also  assumed  that  the  rotor  of 
a  RVR  helicopter  will  operate  in  or  near  to  an  auto-rotative  condition  (small  amount  of 
rotor  power  at  low  tip  path  plane  angle)  with  auxiliary  propulsion  of  the  vehicle.  [Ewans, 
1973]  While  in  this  flight  condition  the  stability  axis  of  the  aircraft  is  moved  away  from 
the  center  of  the  rotor  disk  to  a  location  closer  to  that  of  a  typical  fixed  wing  aircraft. 
This  relieves  the  lateral  and  longitudinal  cyclic  pitch  constraints  since  the  rotor  disk  will 
seek  an  equilibrium  condition  based  primarily  on  collective  pitch.  The  final  condition 
that  did  not  allow  the  RVR  system  to  trim  at  high  speed  was  rotor  drag.  With  the  RVR 
system  in  full  operation  ( //  =  2  and  negative  pitch  on  the  retreating  side  of  the  rotor  disk) 
conventional  helicopter  rotor  drag  analysis  does  not  work.  Fortunately,  having  the  rotor 
in  a  near  auto-rotative  state  allows  a  simple  but  critical  drag  assumption;  over  the  entire 
rotor,  the  integrated  effect  of  the  tilt  of  the  lift  vector  at  each  blade  element  is  sufficient  to 
overcome  the  integrated  effect  of  the  drag  at  all  of  the  blade  elements.  [Prouty,  1986] 
Only  a  small  amount  of  power  is  required  to  overcome  the  drag  and  provide  rotor 
stability.  These  assumptions  are  highly  interdependent.  All  of  the  conditions  must  be 
satisfied  concurrently  whenever  RVR  trim  calculations  are  applied. 

Since  having  the  correct  advance  ratio,  jj. ,  is  critical  to  successful  RVR  operation 
JANRAD  2P_RVR  has  a  built  in  calculator  which  computes  the  optimal  jj.  for  the  design 
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cruise  airspeed  such  that  the  entire  retreating  blade  is  in  reverse  flow.  This  optimal  jj. 
value  ensures  the  reverse  flow  area  around  the  retreating  blade  is  suffieient  to  support 
RVR  aerodynamic  operation  and  RVR  based  cyclic  pitch,  discussed  below. 

3,  RVR  Based  Cyclic  Pitch 

One  flnal  element  had  to  be  added  to  the  RVR  system;  RVR  based  eyelie  piteh. 
The  use  of  HHSC  piteh  is  required  for  advance  ratios  greater  than  1.0.  [Zientek,  2001] 
Using  1/rev  eyelie  piteh  alone,  the  rotor  cannot  achieve  the  required  large  negative  piteh 
in  the  270  degrees  blade  sector,  yet  maintain  positive  pitch  everywhere  else.  See  Figure 
29. 


InHow  Angle  (deg) 


Figure  29.  Inflow  Angle  versus.  Blade  Station.  [From;  Zientek,  2001] 

The  2/rev  RVR  based  eyelie  piteh  can  obtain  the  proper  blade  angle  of  attaek  on 
the  retreating  side  without  the  aft  tilt  normally  assoeiated  with  windmilling  operation. 
This  differs  from  eonventional  HHSC.  The  equation  for  0  takes  the  form 

^rvr  =  ^0  +  ^0  sin(t^)  -  oos(^^)  -  sin(^^)  +  6^^  cos{2y/) 

The  negative  signs  preceding  the  6^^  and  6^^  terms  are  required  to  move  the 

maximum  eyelical  inputs  to  the  retreating  sections  of  the  rotor  disk.  Figure  30  eompares 
conventional  HHSC  based  eyelie  pitch  and  RVR  based  eyelie  pitch. 
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Figure  30.  Conventional  FIFISC  Cyclic  Pitch  (Left)  Compared  to  RVR  Based  Cyclic 

Pitch. 

D.  EXPERIMENTAL  RESULTS 

All  of  these  assumptions  are  enforced  within  JANRAD  2P_RVR  by  implementing 
user  enabled  logic  gates.  The  logic  gates  allow  JANRAD  2P_RVR  to  operate  in  either 
conventional  or  RVR  based  cyclic  pitch.  JANRAD  2P_RVR  only  contains  one  non-user 
driven  restriction  for  RVR  based  cyclic  pitch  employment,  airspeed  must  be  greater  than 
200  knots.  Figures  31  through  35  show  JANRAD  2P_RVR  output. 

Figure  3 1  clearly  shows  the  effects  of  the  RVR  Based  Cyclic  Pitch.  The  negative 
blade  angles  at  rotor  azimuth  positions  between  230  degrees  and  320  degrees  are  crucial 
to  the  development  of  positive  lift.  Also  note  that  collective  pitch  has  been  reduced  to 
zero  through  JANRAD  2P_RVR’s  trim  constraints,  complementing  the  assumption  that 
the  rotor  system  has  entered  a  near  autorotative  state. 
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Figure  3 1 .  Nature  of  RVR  Based  Cyclic  Pitch. 


Figure  32  shows  the  Mach  number  distribution  around  the  rotor  disk.  Note  that 
the  retreating  side  of  the  blade,  approximate  blade  azimuth  positions  of  200  degrees  to 
340  degrees,  is  fully  entrenched  in  reverse  flow.  It  is  the  combination  of  this  reverse 
flow,  the  RVR’s  local  angle  of  attack,  and  Cl  that  allow  for  positive  lift  distribution. 
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Figure  33  shows  Cl  distribution  viewed  as  a  contour  plot.  One  can  immediately 
see  that  the  disk  is  trying  to  achieve  equilibrium.  HHSC  contour  lines  are  replaced  with 
the  RVR  based  cyclic  pitch  contour  lines.  The  RVR  based  cyclic  pitch  builds  upon  the 
principles  of  FUTSC  but  shifts  the  emphasis  on  stall  control  from  the  forward  and  aft 
portions  of  the  disk  to  the  left  and  right  sides  of  the  disk.  The  shift  of  stall  control  to  the 
lateral  sections  of  the  disk  is  critical  to  take  full  advantage  of  the  traits  of  the  RVR.  The 
contour  lines  in  Figure  33  clearly  illustrate  that  the  RVR  based  cyclic  pitch  is  functioning 
as  predicted.  Also  apparent  in  Figure  33  are  the  high  Cl  values  in  the  forward  right 
sector  of  the  disk.  This  phenomenon  is  discussed  in  the  following  paragraphs.  The 
distinct  left  and  right  side  loading  is  further  illustrated  in  Figures  34  and  35. 
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CL  Distribution  At  300.3893  Kts,  advance  ratio=2.0003.  2P  input=9  degs.  TPP  angle=3  degs.  3P  input=0  degs.  3P  ofFset=0  degs 


90 


Figure  33.  Cl  Distribution  Contour  Plot. 


Figure  34  shows  the  radial  distribution  of  blade  lift  at  the  most  important  rotor 
azimuth  locations.  At  the  270  degree  location  positive  lift  is  observed  even  though 
the  blade  is  total  reverse  flow.  As  stated  earlier,  the  positive  air  load  is  a  result  of  the 
RVR’s  unique  aerodynamic  behavior  given  the  proper  operating  conditions. 
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Figure  34.  Blade  Position  versus  Air  Load  at  Critical  Rotor  Azimuth  Locations  (0, 

90,180,  and  270  Degrees). 


Figure  35  is  a  three  dimensional  representation  of  the  lift  distribution  around  the 
disk.  Note  the  predominant  lift  ‘hump’  centered  at  the  disk’s  270  degree  azimuth 
location,  further  proof  the  RVR  based  cyclic  pitch  is  performing  as  predicted.  Also  note 
the  ‘spike’  centered  at  approximately  ^=160  degrees.  This  ‘spike’  is  the  result  of  high 
blade  angle,  and  therefore  high  blade  coefficient  of  lift.  The  high  blade  angle  is  due  to 
the  cumulative  effects  of  the  cyclic  inputs  within  that  azimuth  quadrant.  Although 
undesirable,  this  phenomenon  does  not  effect  the  disk’s  ability  to  trim  within  the 
computational  boundaries  established  in  JANRAD  2P_RVR. 
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Figure  35.  Lift  Distribution  Surface  Plot. 


This  area  of  abnormal  lift  distribution  will  certainly  cause  vibrations  and  high 
blade  stress  as  each  rotor  blade  is  forced  to  pass  through  the  effected  azimuths.  In  fact, 
Ewans  encountered  similar  unbalanced  lift  forces  in  his  tests  [Ewans,  1975]  that  made  it 
impossible  to  complete  his  entire  test  program.  In  an  effort  to  smooth  the  lift  distribution 
a  3P  signal  (4  degrees)  was  added.  The  resulting  lift  distribution  is  shown  in  Eigure  36. 
In  this  figure  three  distinct  humps  are  visible  and  the  overall  lift  distribution  is  generally 
equalized  between  the  three  humps.  Application  of  additional  higher  harmonic  inputs 
would  continue  to  suppress  unwanted  lift  spikes.  Detailed  study  of  rotor  system  behavior 
using  harmonics  higher  than  3P  is  beyond  the  scope  of  this  study  and  would  require 
further  analysis  in  future  research. 
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Lift  Distribution  At  300.3893  Kts,  advance  ratio=2-0003. 2P  input=3  degs,  3P  input=4  degs,  3P  offset=0  degs 


Figure  36.  Lift  Distribution  Surface  Plot  with  3P  Flarmonic  Cyclic  Pitch. 

JANRAD  2P_RVR  has  shown  that  a  RVR  system  with  tailored  2P  cyclic  pitch 
does  trim  under  the  assumptions  used  in  the  study.  The  ability  to  trim  allows  for 
sustained  high  speed  aerodynamic  flight  without  incurring  high  drag  penalties  or  autogiro 
type  behaviors.  The  ability  to  maintain  aerodynamic  flight  with  a  rotor  system  above  200 
knots  allows  designers  to  expand  the  performance  and  usability  of  current  helicopters. 
The  next  section  explores  one  possible  way  to  employ  a  complete  RVR  system  to  fulfill 
an  emerging  vertical  heavy  lift  requirement. 
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IV.  APPLICATION  OF  THE  RVR  AND  ACCOMPANYING 
CONTROL  SYSTEM  TO  THE  JOINT  VERTICAL  HEAVY  LIFT 

(JVHL)  AIRCRAFT 

A,  FIRST  ITERATION  DESIGNS  TO  SATISFY  REQUIREMENTS 

In  an  effort  to  satisfy  operational  requirements  defined  in  the  first  ehapter,  two 
‘new-start’  aircraft,  the  Joint  Heavy  Lift  (JHL)  and  Reverse  Velocity  Rotor  (RVR),  were 
the  subjects  of  simultaneous  preliminary  design  analysis.  For  more  information  on  these 
advanced  designs  the  reader  is  referred  to  the  references  entitled  “Joint  Heavy  Lift 
Expeditionary  Warfare  Aircraft  Solution”  [Wood  and  Aaron,  2003]  and  “Reverse 
Velocity  Rotor  Approach  for  Design  of  a  STOVL  Heavy  Lift  Transport  for  Expeditionary 
Warfare.”  [Wood  and  Van  Riper,  2003]  The  JHE  is  intended  to  be  a  long  range,  heavy 
lift  aircraft  to  support  Marine  Corps,  joint,  and  coalition  force  operations  ashore  up  to  200 
NM  inland  in  a  forcible  entry  environment.  The  JHE  is  a  Joint  Strike  Eighter  (JSE)  lift- 
fan  variant  modified  for  the  increased  load-carrying  requirement.  Earge  lift  fans,  each 
producing  up  to  60,000  pounds  of  vertical  thrust,  are  embedded  in  the  wing  sections  of 
the  aircraft  to  provide  Vertical  Take-off  and  Eanding  (VTOE)  capability  while  thrust 
vectored  turbofan  engines  provide  forward  thrust  while  in  forward  flight.  The  RVR 
aircraft  is  a  compound  helicopter  equipped  with  a  conventional  anti-torque  tail  rotor  and 
large  turboprop  engines  to  provide  the  auxiliary  thrust  required  for  high-speed  forward 
flight.  The  rotor  system  is  an  eight  bladed  (110  foot  diameter),  fully  articulated,  foldable 
system  with  each  blade  incorporating  the  RVR  airfoil  cross  section.  The  rotor  system  is 
driven  by  a  lightweight  variable  speed  transmission  that  allows  RPM  to  be  set  at  either 
100%  or  50%  of  normal  operating  RPM  based  on  flight  mode.  The  anti-torque  system  is 
a  traditional  six  bladed  pusher  tail  rotor  mounted  to  the  tail  vertical  pylon.  Auxiliary 
propulsion  is  provided  by  two  propellers  (fifteen  foot  diameter),  each  driven  by  one  of 
the  two  wing-mounted  turboprop  engines.  Eigures  37  and  38  show  each  design. 
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Figure  37.  RVR  Based  Aircraft. 


Figure  38.  JFIL  Aircraft. 


Both  aircraft  are  intended  to  operate  from  ships  of  the  amphibious  task  force, 
specifically  the  Amphibious  Assault  and  Logistics  Support  variants  of  the  family  of 
ExWar  platforms  currently  under  design  by  the  Total  Ship  Systems  Engineering  (TSSE) 
curriculum  at  the  Naval  Postgraduate  School.  The  aircraft  are  also  designed  to  be 
compatible  with  legacy  force  ships,  i.e.  EHD,  EHA  (R),  and  MPF  (E)  class  ships  to  the 
maximum  extent  possible. 

As  might  be  expected  these  evolutionary  designs  became  the  subject  of  much 
debate  and  criticism.  Concerns  about  the  RVR  ranged  from  the  expected  size,  weight, 
and  complexity  of  a  variable  speed  main  transmission  able  to  handle  an  aircraft  with  a 
take-off  gross  weight  of  over  100,000  pounds  to  the  employment  of  a  tail  rotor  almost  as 
large  as  a  UFI-60  main  rotor  blade.  The  JFIE  suffered  from  lift  fan  and  wing  integration 
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challenges  and  a  downwash  veloeity,  the  vertieal  wind  produeed  as  the  aircraft  hovers, 
whieh  eould  praetieally  lift  eonerete  bloeks  off  the  ground  and  send  them  flying  several 
yards  away.  Clearly,  a  more  amiable  solution  had  to  exist. 

B,  THE  JOINT  VERTICAL  HEAVY  LIFT  AIRCRAFT 

Ongoing  research  has  produced  a  new  eoneeptual  design  that  eombines  the  best 
traits  of  the  RVR  and  JHL  aircraft  while  greatly  redueing  the  negative  aspeets  of  eaeh 
design.  Although  the  JHL  did  not  employ  any  RVR  or  HHSC  teehnologies  that  form  the 
eore  of  this  study,  it  provides  the  framework  for  the  FI  19  (JSF)  Shaft  Driven  Lift  Fan  and 
Propulsion  module  integration.  This  teehnology  provides  the  key  to  overeoming  the  most 
ehallenging  design  point  of  the  aircraft;  sustained  hover  and  low  speed  flight  in  harsh 
environments.  The  proposed  hybrid  design,  labeled  the  Joint  Vertieal  Heavy  Lift  (JVHL) 
uses  a  RVR  based  rotor  system  and  drive  train  with  three  F-35  V/STOL  type  Shaft 
Driven  Lift  Fan  and  Propulsion  Modules.  A  direct  jet  anti-torque  thruster  is  used  in  lieu 
of  a  eonventional  anti-torque  tail  rotor.  JANRAD  2P_RVR  was  used  exclusively  for  the 
JVHL’s  rotor  performance  caleulations.  A  detailed  description  of  JANRAD  2P_RVR 
and  user  instruetions  are  provided  in  Chapter  V.  The  proposed  eonfiguration  of  the 
JVHL  is  shown  in  Figures  39-42.  Dimensions  shown  in  Figure  39  are  in  feet. 


Figure  39.  JVHL  Three  View. 
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^  C-130J  sized  fuselage 


Figure  40.  Top  View  of  JVHL. 
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Figure  41. 


Perspective  View  of  JVFIL. 


Figure  42  shows  the  JVFIL  in  hover  configuration  with  60,000  lbs  of  lift  from  the 
main  rotor,  18,000  pounds  from  each  fuselage  mounted  FI  19  module  and  18,000  pounds 
from  the  tail  mounted  shaft  driven  lift  fan,  for  a  total  of  1 14,000  lbs  of  vertical  lift. 
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Figure  42.  JVFIL  in  Flover  Configuration. 

Immediate  benefits  of  using  the  FI  19  propulsion  modules  are  elearly  evident. 
Most  important,  the  propulsion  module  offsets  the  hover  and  low  speed  performance  of 
the  RVR  system  by  providing  a  great  deal  of  vertical  thrust.  Current  FI  19  engine 
configurations  would  provide  ample  power  for  all  flight  modes  of  the  JVHL.  In  fact,  the 
three  FI  19  modules  would  provide  such  a  surplus  of  power  that  each  module  could  be 
down-rated  by  up  to  30%.  This  modification  would  undoubtedly  reduce  fuel 
requirements  and,  more  importantly,  reduce  engine  fatigue. 

1,  Fuselage  Modules 

The  fuselage  mounted  modules  are  installed  so  that  the  drive  shafts  leading  to  the 
main  transmission  have  a  straight  path  and  direct  jet  exhausts  do  not  have  to  be  ducted 
around  the  fuselage.  When  in  forward  flight  the  two  fuselage  mounted  modules  provide 
ample  thrust  so  the  JVHL  can  achieve  the  design  speed  of  300  knots  while  at  an  altitude 
of  10,000  feet  MSL.  The  large  compound  wing  no  longer  has  to  accommodate  turboprop 
engines  or  complicated  shafting.  The  wing  can  now  have  a  full  complement  of  control 
surfaces  and  even  be  modified  to  carry  extra  fuel.  Of  course  the  wing  will  still  have  a 
negative  effect  on  hover  performance  but  the  added  vertical  thrust  provided  by  the  direct 
jet  thrusters  and  the  aft  lift  fan  help  reduce  the  wing’s  adverse  effects. 

2.  Tail  Module 

The  tail  module  eliminates  all  tail  rotor  components  and  replaces  them  with  a  self- 
contained  highly  reliable  vertical  lift  augmentation  and  anti-torque  solution  while 
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significantly  reducing  drag  in  high-speed  flight.  The  vectored  thrust  nozzle  is  collective 
and  pedal  coupled  to  ensure  maximum  anti-torque  thrust  and  directional  control  is 
delivered  proportional  to  inputs.  While  hovering,  the  module  would  also  provide 
approximately  20,000  lbs  vertical  lift  via  the  built  in  Lift  Fan.  This  lift  would  help  offset 
main  rotor  lift  and  power  requirements  and  extend  the  CG  range  of  the  aircraft. 

3,  Preliminary  Design  Data 

Another  benefit  of  the  main  rotor,  direct  jet  thruster,  and  left  fan  configuration  is 
that  the  downwash  will  be  less  than  half  that  of  the  JHL  and  the  direct  jet  thruster  output 
will  be  almost  immediately  cooled  by  the  entrained  flow  circulating  through  the  main 
rotor  system.  One  of  the  most  important  advantages  of  using  the  FI  19  module  is  that 
now  the  JVHL  will  have  a  measurable  power  and  performance  margin  that  will  enable 
the  aircraft  to  hover  with  a  full  load  on  a  Fligh/FIot  day  (4000  ft  MSL  and  95  degrees  F). 
Critical  design  data  are  shown  in  Tables  1  and  2. 


High  Hot.  4000ft.  95degree  F,  rho 
is  0.00192 

Standard  Dai 

Nointal  drive 

Normal  drive 

Parameters 

Notation 

JVHL  (hover) 

JVHL  (fwd 
flight) 

JVHL  (hover) 

JVHL  (fwd 
flight) 

Weiqhts  (lbs) 

Empty 

W, 

62865 

62865 

62875 

62875 

Payload 

37500 

37500 

37500 

37500 

Fuel  Capacity 

W<.,i 

20000 

20000 

20000 

20000 

Empty/Gross  ratio 

W./GW 

0.52 

0.52 

0.52 

0.52 

Gross  weight 

GW 

120365 

120365 

120375 

120375 

Rotor/compound  winq  proportion 

Lift  from  direct  jet/aft 
fan 

LF 

50000 

50000 

%  of  lift  for  rotor 

100% 

20% 

100% 

20% 

of  lift  for  winq 

0% 

80% 

0% 

80% 

Rotor  parameters 

Thrust  (lbs) 

T  =  %L„i„*GW 

70365 

'24073 

70375 

24075 

Radius  (ft) 

R 

40 

40 

40 

40 

Chord  (ft) 

c 

3 

3 

3 

3 

Ho.  of  blades 

b 

8 

8 

8 

8 

FWD  Speed  (kts), 

Okts  means  hover 

V.-j 

0 

300 

0 

300 

Omeqa  (RPM) 

Omeqa 

150 

60 

150 

60  ' 

Figure  of  merit 

FM 

0.8 

0.8 

0.8 

0.8 

Disc  area  (ft''2) 

A  =  pi*R''2 

5027.20 

5027.20 

5027.20 

5027.20 

Solidity 

Sigma  = 
(b*c)/(pi*R) 

0.191 

0.191 

0.191 

0.191 

Table  1.  Basic  JVFIL  Design  Data. 
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Hi9k  Hot,  4000ft,  95de9rec  F,  rko  is 
0.00132 

Steederd  Def 

Nora^l  drive 

Nornel  drive 

Parameters 

Notation 

JVHL  (hover) 

JVHL  (fwd 
flight) 

JVHL  (hover) 

JVHL  (fwd 
flight) 

Tip  speed  at  0  deg 
rotor  azimuth 

Vllr.i  = 

Omeqa*R 

628.40 

251.36 

628.40 

251.36 

Tip  speed  at  90  deg 
rotor  azimuth 

V  lip, 11  = 

Omeqa'R  •  Vf-i 

628.40 

758.36 

628.40 

758.36 

Tip  speed  at  180  deg 
rotor  azimuth 

V  = 

Omeqa'R 

628.40 

251.36 

628.40 

251.36 

Tip  speed  at  270  deg 
rotor  azimuth 

Vi:p.z7i  = 

Omeqa'R  -  V«.j 

628.40 

-255.64 

628.40 

-255.64 

Mai  Mach  lip 

Mai  Mi:, 

0.57 

0.68 

0.57 

0.68 

Disc  loading  fpsfl 

DL  =  T/ A 

14.00 

4.79 

14.00 

4.79 

Induced  velociti  f^ps) 

V:  = 

sgrtrOL/2'rho1 

60.37 

Not  applicable 

54.25 

Not  applicable 

Advance  ratio 

mu  = 

V«.j  fOmeqa'R 

0.00 

w 

2.01 

0.00 

w 

2.01 

P:.j...jrHP| 

P:.j...j(HP)  = 
T’V:  1  550 

7724.02 

Not  applicable 

6941.94 

Not  applicable 

1st  approz: 

P,.i..ifHPl 

. . (HP)  = 

P:.j...j  f  FM 

9655.03 

Not  applicable 

8677.43 

Not  applicable 

Thrust  coed.  1  sigma 

Ct  = 

T 1  A'tho*(omega 
'R1*2'sigma 

0.088 

Not  applicable 

0.078 

Not  applicable 

Torgue  coed.  1  sigma 

C*  1  sigma 
(from  Proutlg 
teit  page  90-92) 

0.0270 

Not  applicable 

0.020 

Not  applicable 

Torgue 

Q  = 

786083 

Not  applicable 

721183 

Not  applicable 

2nd  approi: 

P-i„ifHPl 

calculated  from 

P  = 

Q*omegaf550 
where  Q  is  from 
Cqfsigma  curve 
in  teit  page  90- 
92 

22450 

Not  applicable 

20597 

Not  applicable 

plus  33X 

Final  P  reguired 

29859 

Not  applicable 

27394 

Not  applicable 

Table  2.  JVHL  Rotor  and  Power  Calculations. 

C.  RVR  ANALYSIS 

The  JVHL  uses  two  emerging  technologies,  the  RVR  system  and  the  FI  19 
propulsion  module.  All  other  major  components  of  the  aircraft  are  based  on  fairly  mature 
technologies  and  processes.  Compound  helicopter  design  is  well  understood,  the  CH- 
53E  rated  transmission  and  drive  line  technology  needed  to  satisfy  main  rotor  lift 
requirements  has  been  refined  by  the  helicopter  industry,  and  use  of  composites  in 
fuselage  design  practices  is  becoming  more  commonplace  with  each  passing  year.  The 
only  remaining  design  obstacle  is  the  RVR  rotor  system. 

A  detailed  aeroelastic  and  dynamic  analysis  of  the  RVR  system  is  beyond  the 
scope  of  this  study  but  a  preliminary  aerodynamic  analysis  was  conducted  to  ensure  the 
RVR  system  would  be  viable  at  high  airspeeds.  JANRAD  version  lP2P3Prvr  was  used 
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to  check  rotor  system  performance.  As  discussed  earlier,  correctly  manipulating  the 
advance  ratio  is  critical  since  the  RVR  system  relies  on  the  entire  retreating  side  of  the 
rotor  system  being  fully  entrenched  in  reverse  airflow.  For  the  JVHL  an  advance  ratio  of 
2  is  used  at  the  aircraft’s  cruise  speed  of  300  knots.  To  achieve  this  advance  ratio  the 
rotational  speed  of  the  rotor  system  is  reduced  from  150  rpm  to  60  rpm.  This  high 
advance  ratio  causes  sufficient  reverse  airflow  on  the  retreating  side  of  the  rotor  disk, 
trailing  edge  to  leading  edge,  to  produce  lift.  The  RVR  is  coupled  with  a  RVR  unique 
two  per  revolution  (2/rev)  cyclic  pitch  signal.  See  Figures  43-45. 


Rotor  Azimutii  Angle  v.  RVR  Harmonic  Cyclic  Pitch  for  300  3893  Kts.  advance  ratlo=20003.  2P  input=9  degs 


Rotor  Azimuth  Angle  (degs) 


Figure  43.  RVR  Based  Cyclic  Pitch  for  the  JVHL. 
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180 

Mach  Number  Distribution  At  300.3893  Kis,  advance  ratio=2.0002 


Figure  44.  Mach  Number  Distribution  at  Cruise  Speed  (300  kts), Retreating  Side  of 

Blade  is  in  Reverse  Flow. 


Figure  45.  Cl  Distribution  at  Cruise  Speed  (300  knots)  for  the  JVHL. 

Since  the  large  compound  wing  has  unloaded  the  rotor  (at  cruise  speed  the  wing  is 

carrying  80%  of  aircraft  weight  and  the  rotor  is  carrying  the  remaining  20%),  the  rotor 
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disk  can  operate  in  an  autorotative  state  with  small  cyclic  input  and  no  wind  milling. 
Without  this  wind  milling  behavior  there  is  almost  no  drag  penalty.  These  figures 
illustrate  the  behavior  of  this  rotor  system  is  directly  in  agreement  with  predicted  RVR 
behavior.  The  aircraft  is  able  to  operate  at  high  cruise  speed  without  the  typical  negative 
effects  of  the  rotor  disk. 
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V.  JANRAD  USER  INSTRUCTIONS 


A.  GENERAL 

Program  installation,  guidance  for  input  requirements,  and  rotor  performanee 
output  are  diseussed  in  the  subsequent  seetions.  These  instruetions  assume  that  the  reader 
has  a  working  knowledge  of  MATLAB  release  12  and  higher. 

Joint  Army  Navy  Rotoreraft  Analysis  and  Design  (JANRAD)  version  lP2P3Prvr 
was  developed  to  allow  introduetion  of  2/rev  eyelie  inputs,  3/rev  eyelie  inputs,  eompound 
wing  sizing,  uneonventional  vertieal  lift  augmentation,  and  employment  of  the  Reverse 
Veloeity  Rotor  (RVR)  system  and  subsequent  RVR  speeifie  2P  eontrol  inputs.  This 
version  also  ineludes  a  eomprehensive  plotting  eapability  that  gives  the  user  high  fidelity 
2-D,  3-D,  and  eontour  plots  to  aid  in  rotor  performanee  analysis. 

JANRAD  version  lP2P3Prvr  was  built  using  the  original  souree  eode  found  in 
JANRAD  version  1  [Nieholson,  1993].  This  approaeh  ensured  the  original  trim  routine 
that  forms  the  heart  of  program  was  unaltered  by  subsequent  undoeumented  updates. 

The  program  was  developed  and  tested  on  a  Pentium  4  based  laptop  running 
MATLAB  release  13  under  Windows  XP  Professional.  The  laptop  was  equipped  with 
512  megabytes  of  RAM  and  a  twenty  gigabit  fixed  drive. 

B,  INSTALLATION 

To  install  JANRAD  2P_RVR,  ereate  a  subdireetory  named  JANRAD  and  install 
the  following  JANRAD  2P_RVR  modules:  janrad.m,  trim.m,  dmeale.m,  hh02eled.m, 
ool2eled.m,  perfm,  plotter.m,  polarl.m,  rvreled.m,  sel095r8eled.m,  threale.m,  tmeale.m, 
and  vrl2eled.m.  Listings  of  the  modules  ean  be  found  in  subsequent  appendiees. 

Onee  this  installation  is  eomplete,  open  MATLAB  and  type  JANRAD  (with  the 
proper  drive  designation)  in  the  Current  Direetory  dialogue  box.  All  of  the  modules  must 
be  installed  in  the  seleet  direetory  on  the  eomputer’s  fixed  drive  for  the  program  to 
funetion.  The  program  will  not  work  if  the  user  attempts  to  run  the  modules  from  a 
remote  drive,  eompaet  disk,  or  network.  The  program  is  now  ready  to  run. 
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C.  PROGRAM  EXECUTION 
1.  Input 

To  run  JANRAD  2P_RVR  beta,  at  the  MATLAB  command  prompt,  type 
‘janrad’.  At  the  Title  page  enter  a  1  if  you  want  to  perform  analysis  on  a  conventional 
helicopter  or  2  if  you  want  to  perform  analysis  on  a  compound  helicopter.  RVR  based 
analysis  is  only  allowed  if  you  select  the  compound  helicopter  option.  See  Figure  46. 


Figure  46.  Title  Page. 

At  the  File  Selection  page,  enter  a  1  if  you  want  to  edit  a  previously  stored  data 
file,  otherwise  enter  a  2.  If  you  chose  to  edit  a  file,  you  will  be  prompted  on  the  next 
page  to  enter  the  file  name  without  the  extension.  If  JANRAD  2P_RVR  is  unable  to  find 
the  specified  file  you  will  be  prompted  to  try  again.  Ensure  your  spelling  is  correct  and 
you  are  in  the  proper  directory.  Once  the  file  is  loaded  the  edit  menu,  as  seen  in  Figure 
47,  will  be  displayed.  Check  to  make  sure  the  menu’s  title  corresponds  to  the  type  of 
helicopter  you  are  analyzing,  either  conventional  or  compound. 
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COHPOUMD  HELICOPTER 

EDIT  MENU 

1. 

pressure  altitude 

2. 

temperature 

3. 

airspeed 

4. 

gross  weight 

5. 

numlier  of  blades 

6. 

blade  radius 

7. 

blade  chord 

8. 

hinge  offset 

9. 

blade  grip  length 

10. 

blade  tTnist 

11. 

blade  wieght 

12. 

number  of  blade  elements 

13. 

rotational  velocity 

14. 

if  azimuth  sectors 

15. 

lift  curve  slope 

IS. 

airfoil 

17. 

collective  pitch 

10. 

flatplate  area 

19. 

vert  projected  area 

20. 

horizontal  tail  area 

21. 

horizontal  tail  span 

22. 

horizontal  tail  CL 

23. 

horizontal  tail  CDo 

24. 

vertical  tail  area 

25. 

vertical  tail  span 

26. 

vertical  tail  CL 

27. 

vertical  tail  CDo 

28. 

2P  input 

29. 

3P  input 

30. 

auxiliary  thrust 

31. 

Tiling  area 

32. 

wing  span 

33. 

Tiring  CL 

34. 

wing  CDo 

35. 

Tiring  efficiency  factor 

0. 

NO  CHANGES 

Select  the  patameter  to  change 

=  1 

Figure  47.  Compound  Helicopter  Edit  Menu  Page. 


Select  the  number  of  the  parameter  you  wish  to  change.  The  current  value  for  that 
parameter  will  be  displayed  followed  by  a  prompt  to  enter  the  new  value.  If  you  chose 
not  to  change  the  value,  press  <ENTER>  and  the  previous  value  will  be  maintained. 
Once  a  new  value  is  entered  or  <ENTER>  is  pressed  the  edit  menu  will  be  displayed 
again.  Choose  the  number  of  the  next  parameter  to  change  or  0  to  exit  the  programs 
direct  edit  mode.  If  you  chose  to  create  a  new  data  file,  prompts  for  individual 
parameters  will  be  displayed.  Enter  the  value  of  each  parameter  when  prompted. 

Enter  the  values  for  the  specified  parameters  regardless  of  editing  or  creating  a 
data  file  in  the  following  manner; 

1 .  PA  -  Enter  the  pressure  altitude  in  feet 

2.  temperature  -  Enter  the  temperature  in  degrees  Eahrenheit 

3.  airspeed  -  Enter  forward  airspeed  in  knots.  JANRAD  automatically 
converts  airspeed  to  feet/second. 

4.  GW  -  Enter  aircraft  gross  weight  in  pounds. 

5.  number  of  blades  -  Enter  the  number  of  blades  used  in  the  rotor  system. 
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6.  radius  -  Enter  the  rotor  blade  radius  in  feet.  Measure  the  radius  from  the 
center  of  the  rotor  hub  to  the  tip  of  the  bade. 

7.  chord  -  Enter  the  average  rotor  blade  chord  in  feet. 

8.  hinge  offset  -  Enter  the  effective  hinge  offset  in  feet. 

9.  blade  grip  -  Enter  the  distance  from  the  center  of  the  rotor  hub  to  the 
beginning  of  the  aerodynamic  portion  of  the  rotor  blade  in  feet. 

10.  blade  twist  -  Enter  rotor  blade  twist  in  degrees.  You  can  enter  blade  twist 
as  either  a  positive  or  negative  value.  JANRAD  uses  the  absolute  value  of 
twist. 

1 1 .  blade  weight  -  Enter  the  weight  of  a  single  rotor  blade  in  pounds. 

12.  rotor  blade  elements  -  Enter  the  number  of  rotor  blade  radial  elements. 
Multiples  of  20  are  recommended,  i.e.  20,  40,  80.  As  the  number  of  blade 
elements  increases  program  execution  is  slowed.  JANRAD  adds  one 
additional  element  to  account  for  tip  loss. 

13.  rotational  velocity  -  Enter  the  rotor  rotational  velocity  in  radians/second. 

14.  azimuth  sectors  -  Enter  the  number  of  rotor  disk  azimuth  sectors.  This 
number  must  match  the  number  of  rotor  blade  elements  if  JANRAD 
created  plots  are  to  be  used.  The  number  may  be  different  if  no  plots  are 
required. 

15.  lift  curve  slope  -  Enter  the  average  slope  of  the  linear  portion  of  the  rotor 
blade’s  airfoil’s  lift  curve,  where  the  lift  curve  is  Cl  versus  alpha. 

16.  rotor  blade  airfoil  -  Eour  airfoils  are  given  for  use  with  conventional 
helicopters  and  five  airfoils  are  given  for  use  with  compound  helicopters. 
Enter  the  number  that  corresponds  to  your  desired  airfoil. 

17.  collective  pitch  -  Enter  the  estimated  collective  pitch  at  0.7  r/R  in  degrees. 
JANRAD  automatically  converts  the  vale  to  radians. 

18.  flat  plate  area  -  Enter  the  aircraft  equivalent  flat  plate  area  in  square  feet. 
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19.  vertical  projected  area  -  Enter  the  vertical  projected  area  of  the  aircraft  in 
square  feet. 

20.  horizontal  tail  area  -  Enter  the  horizontal  tail  area  in  square  feet.  If  no 
horizontal  tail,  enter  a  0. 

2 1 .  horizontal  tail  span  -  Enter  the  horizontal  tail  span  in  feet. 

22.  horizontal  tail  Cl  -  Enter  the  expected  lift  coefficient  for  the  horizontal 
tail. 

23.  horizontal  tail  Cdo  -  Enter  the  horizontal  tail  profile  drag  coefficient. 

24.  vertical  tail  area  -  Enter  the  area  of  the  vertical  tail  in  square  feet.  If  no 
vertical  tail,  enter  a  0. 

25.  vertical  tail  span  -  Enter  the  span  of  the  vertical  tail  in  feet. 

26.  vertical  tail  Cl  -  Enter  the  expected  lift  coefficient  for  the  vertical  tail. 

27.  vertical  tail  Cdo  -  Enter  the  vertical  tail  profile  drag  coefficient. 

28.  2P  longitudinal  input  -  Enter  the  2/revolution  fixed  cyclic  input. 

29.  3P  input  -  Enter  the  3/revolution  fixed  cyclic  input.  Enter  0  for  no  3P 
input.  Regardless  of  3P  input  you  must  enter  an  offset  value  between  0 
and  119  degrees. 

30.  auxiliary  thrust  -  Enter  the  expected  auxiliary  thrust  in  pounds. 

3 1 .  wing  area  -  Compound  helicopter  edit  menu  only.  Enter  the  wing  area  in 
square  feet. 

32.  wing  span  -  Compound  helicopter  edit  menu  only.  Enter  the  wing  span  in 
feet.  If  you  included  the  segment  of  fuselage  encompassed  by  the  wing  in 
the  wing  area  value  include  it  in  the  span  value  also. 

33.  wing  Cl  -  Compound  helicopter  edit  menu  only.  Enter  the  expected  lift 
coefficient  for  the  wing. 

34.  wing  Cdo  -  Compound  helicopter  edit  menu  only.  Enter  the  wing  profile 
drag  coefficient. 
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35.  wing  efficiency  factor  -  Compound  helicopter  edit  menu  only.  Enter  the 
wing’s  expected  efficiency  factor. 

If  you  have  chosen  Compound  Helicopter  Analysis  you  will  be  prompted  to 
answer  three  additional  questions;  1.  Do  you  want  to  fix  the  Tip  Path  Plane?  2.  Do  you 
want  JANRAD  to  set  thrust  equal  to  total  drag?  3.  Do  you  want  JANRAD  2P_RVRto 
use  RVR  based  2/rev  Cyclic  Pitch?  It  is  recommended  that  the  tip  path  plane  be  fixed  for 
compound  helicopters.  If  you  chose  to  fix  the  angle  you  will  be  prompted  to  enter  the 
angle  in  degrees.  If  you  chose  not  fix  the  angle,  the  program  will  modulate  the  angle  to 
achieve  desired  performance  characteristics.  If  you  chose  not  to  set  auxiliary  thrust  to 
total  drag  JANRAD  use  whatever  value  was  entered  in  the  edit  menu.  Only  chose  to  use 
2/rev  cyclic  pitch  if  you  are  performing  high  speed  (greater  than  190  kts)  slowed  rotor 
analysis.  If  Conventional  Helicopter  Analysis  was  selected,  these  questions  are  bypassed 
and  you  are  taken  directly  to  the  save  instructions  page. 

After  all  the  required  parameters  are  entered,  you  will  be  prompted  to  save  the  data 
under  a  user  designated  file  name.  The  file  name  may  be  up  to  twenty  characters  long. 
JANRAD  2P_RVR  will  save  the  tile  to  the  current  directory  as  the  file  name  with  a 
“.mat”  extension.  If  you  were  editing  a  data  tile,  press  <ENTER>  if  you  want  to  save  the 
file  with  the  original  file  name.  Eigure  48  shows  the  Save  Instructions  page. 


***  SAVE  IH3TRUCTI0H3  *** 

A.  Save  the  new  data  to  a  specif iced  file  name. 

B.  Dg  not  use  an  extension  or  quotations. 

C.  Use  the  letter /number  combinations  of  20  characters  or  less. 

D.  The  file  will  be  saved  with  a  ".mat"  extension. 

ex:  dsgn_2 

E.  If  you  made  no  changes,  press  <  Enter  >,  the  file  will 
be  saved  with  the  original  name. 

save  file  as:  | 

Eigure  48.  Save  Instructions  Page. 
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The  Execution  Menu  page  is  displayed  next.  Enter  1  if  you  want  to  proceed  to 
next  step  in  Rotor  Performance  Analysis,  a  2  if  you  would  like  to  go  back  and  edit  your 
data  further  or  if  you  want  to  create  a  new  file.  To  quit  the  program,  select  3. 

After  selecting  Rotor  Performance  Analysis,  JANRAD  2P_RVR  will  ask  you 
several  more  questions  based  on  user  defined  aircraft  design  inputs.  Each  question  is 
discussed  below. 

If  a  Compound  Helicopter  Analysis  was  selected,  JANRAD  2P_RVR  will  display 
the  percent  of  aircraft  weight  carried  by  the  wing.  You  may  choose  to  alter  the 
percentage  of  weight  carried  by  the  wing  or  leave  it  at  the  original  value.  If  you  decide  to 
change  the  value  the  program  will  calculate  and  display  the  new  required  wing  area. 

Regardless  of  Conventional  or  Compound  analysis  JANRAD  will  check  to 
determine  if  hover/low  speed  (airspeed  less  than  1 0  kts)  performance  analysis  is  about  to 
be  performed.  If  hover/low  speed  analysis  has  been  selected  you  will  be  asked  if  there 
are  any  auxiliary  vertical  thrusters,  i.e.  direct  jet  thrusters  or  lift  fans,  installed  on  your 
aircraft.  If  installed,  enter  the  amount  of  thrust  in  pounds  provided  by  these  devices. 

If  a  Compound  Helicopter  Analysis  was  selected  and  you  chose  not  to  set 
auxiliary  thrust  to  total  drag,  JANRAD  2P_RVR  will  check  to  determine  if  a  thrust  deficit 
for  you  desired  flight  condition  exist.  If  a  deficit  exists,  it  will  be  displayed  and  you  will 
be  given  another  chance  to  set  auxiliary  thrust  equal  to  total  drag.  If  you  decide  to 
continue  with  the  original  auxiliary  thrust  value,  a  warning  message  is  displayed  and  the 
program  proceeds  to  the  next  analytical  step. 

If  you  previously  chose  to  employ  2/rev  cyclic  pitch,  JANRAD  2P_RVR  will 
display  the  current  value  of  the  advance  ratio,  p. ,  and  if  required,  suggest  altering  p  to  the 
optimum  value  for  RVR  operation.  Advance  ratio  must  be  greater  than  or  equal  to  1.0  to 
enable  RVR  based  cyclic  pitch.  If  you  allow  the  program  to  change  p  the  rotational 
velocity  of  your  rotor  will  be  changed  to  reflect  the  optimized  p . 
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After  all  these  inputs  have  been  satisfied,  the  program  will  display  an  ‘ANALYSIS 
IN  PROGRESS’  message  and  several  status  messages.  If  your  inputs  fail  to  produced  a 
converged  solution  using  either  conventional  or  RVR  trim  criteria  the  message  shown  in 
Figure  49  will  be  displayed. 

This  configuration  will  not  trim 
Try  a  lower  airspeed  or  a  new  design 

Program  execution  terminated 

???  ***  EWD  OF  PROGRAM  *** 

» _ 

Figure  49.  Program  Termination  Page. 


From  this  page  you  must  restart  JANRAD  2P_RVR  from  the  command  prompt. 
A  ‘ROTOR  TRIMMED’  message  will  be  displayed  when  the  analysis  is  complete. 

2,  Output 

Upon  completion  of  the  analysis,  JANRAD  2P_RVR  will  ask  if  you  want  the 
performance  results  displayed  on  screen.  If  you  answer  ‘yes’  two  screens  displaying 
input  data  and  two  screens  displaying  output  data  will  be  displayed.  The  final  output 
screen  is  displayed  in  Figure  50. 


CALCULATED  DATA 

COHTINUED 

uh60a3 

solidity  (sigma) 

0.082 

Disk  Loading 

0.00 

lbs/ft'"2 

Figure  of  Merit  = 

0.00 

CT/sigma 

0.081 

CQ/sigma 

0.0079 

CH/sigma 

0.0000 

Tip  Mach  numher  of  adv  blade= 

0.808 

Advance  ratio  = 

0.285 

Rotor 

thrust  reguired  (TPP)= 

15366 

lbs 

Rotor  power  reguired= 

1984 

h.p. 

Rotor  torgue= 

41243 

ft- lbs 

press 

any  key  to  continue. . . 

F igure  50.  F ourth  Output  Page . 
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Each  of  these  sereens  is  displayed  one  at  a  time.  As  soon  as  you  advanee  to  the 
next  sereen,  you  eannot  return  to  the  previous  sereen.  A  copy  of  the  four  sereens  is  saved 
to  the  current  directory  with  the  file  name  you  entered  at  the  Save  Instructions  page  and 
eoneatenated  with  a  ‘.prf  extension.  If  you  answer  ‘no’  the  four  screens  will  still  be 
saved  and  you  will  be  asked  if  you  would  like  to  view  any  of  the  JANRAD  2P_RVR 
ereated  plots.  If  you  answer  ‘yes’  the  Available  Plots  page,  seen  below  in  Figure  51,  will 
be  displayed. 


****  Available  Plots  **** 

1.  Rotor  Asimuth  Angle  v.  Cyclic  Pitch 

2.  Blade  AEimuth  Position  v.  Airload 

3.  Rotor  Airload  Distribution 

4.  Blade  Azimuth  Position  v.  Lift 

5.  Rotor  Lift  Distribution 

6.  Angle  of  Attack  Distribution 

7.  Rotor  Tip  Path  Plane  Stall  pattern  and  Thrust  Distribution 

8.  Rotor  Trim  Lissajous  Figure 

9.  Lift  Distribution  (Contour  Plot] 

10.  CL  Distribution  (Contour  Plot) 

11.  Hach  Humber  Distribution  (Contour  Plot) 

Please  choose  a  plot  by  entering  1-11. | 


Figure  5 1 .  Available  Plots  Page. 


Seleeted  plots  will  be  displayed  in  new  windows  labeled  with  the  appropriate 
figure  numbers.  These  plots  may  be  simply  viewed  on  the  screen  or  printed  by  choosing 
the  appropriate  MATFAB  menu  bar  commands.  If  they  are  viewed  on  sereen  it  is 
recommended  that  you  maximize  the  window  for  best  elarity  and  plot  resolution.  When 
printing,  set  your  printer  to  landseape  mode  for  the  best  results.  After  you  have 
eompleted  viewing  the  seleeted  figure,  either  close  or  minimize  the  window.  A  message 
will  be  in  the  eommand  window  asking  if  you  want  to  view  another  plot.  Answer  ‘yes’ 
or  ‘no’.  If  you  chose  ‘yes’  simply  repeat  the  proeess  until  you  have  viewed  all  required 
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plots.  If  you  choose  to  view  plots  2,  4,  or  7,  you  must  manually  position  the  plot  titles. 
To  do  this,  simply  use  the  mouse  to  position  the  cursor  at  the  point  where  you  would  like 
the  title  to  begin  and  click  the  left  button.  Examples  of  each  plot  in  Figures  52-62. 


Figure  52.  Cyclic  Pitch  versus  Rotor  Azimuth  Angle. 
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Figure  53.  Blade  Position  versus  Air  Load. 


Fte  Edit  View  Insert  Tools  Window  Help 


Figure  No.  3 


|D  cS  H  «  t  A  ^  ^  o 

Airload  Distribution  At  120.1557  Ws,  advance  ratio=0.28189, 2P  inpul=4  degs,  3P  input=0  degs,  3P  offset=0  degs 


Figure  54.  Air  Load  Distribution. 


53 


Figure  55.  Blade  Position  versus  Lift. 


File  Edt  View  Insert  Tools  Window  Help 


Figure  No.  5 
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iDciHS  K  A  /<  y  ^  ^  O 

Lifl  Distribution  At  120.1557  Kts,  advance  ratio=0.281B9, 2P  input=4  degs,  3P  input=0  degs.  3P  offset=0  degs 


Figure  56.  Lift  Distribution. 
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Figure  57.  Angle  of  Attack  Distribution. 


Figure  58.  Actual  Versus  Ideal  Rotor  Tip  Path  Plane  Angle  and  Thrust  Distribution. 
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Figure  59.  Lissajou  Figure. 


Figure  60.  Air  Load  Contour  Plot. 
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1  Figure  No.  10 

SEfEil 

RIe  Edit  View  Insert  Tools  Window  Help 
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180 

CL  Distnbution  At  120.1557  Kts,  advance  ratio=0-28189.  2P  input=4  degs.  TPP  angle=534.4577  degs.  3P  input=0  degs.  3P  off5et=0  degs 


Figure  61.  Cl  Distribution  Contour  Plot. 


Figure  62.  Maeh  Number  Distribution. 
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If  you  chose  ‘no’,  the  Performance  Output  Data  Instructions  page  will  be 
displayed.  See  Figure  63. 

***  PERFCiPHAlICE  OUTPUT  DATA  IHSTRUCTI0H3  *** 

A.  Single  value  data  saved  to  a  file  named: 

"uUeOaU.pEf" 

B.  This  is  a  text  file,  use  the  tyoe  command  to  view  the 
file  OE  use  a  text  editoE  to  view/ptint  the  file. 

A  duplicate  file  name  exists,  oe  the  file 
cannot  be  found. 

C.  HatEix  and  vectoE  data  saved  to  a  file  named: 

"uh60a3_p . mat" 

D.  This  is  a  ".mat"  binaty  file,  use  the  load  command 
to  EetEieve  the  data  foE  plotting. 

***  EHD  OF  PERFORHAUCE  *** 

ptess  any  key  to  continue... 

Figure  63.  Performance  Output  Data  Instruetions  Page. 

3,  Saved  Data 

The  veetor  and  matrix  data  generated  by  the  rotor  performanee  module  of 
JANRAD  2P_RVR  is  saved  to  the  eurrent  direetory  with  the  file  name  you  entered  at  the 
Save  Instruetion  page  and  eoneatenated  with  a  ‘_p’  and  a  ‘.mat’  extension.  The 
following  veetors  and  matriees  are  saved. 

1.  r  -  a  vector  containing  the  radii  (in  feet)  to  the  eenter  of  eaeh  blade 
element.  This  is  a  row  veetor  eontaining  n  elements,  where  n  equals  the  number  of  blade 
elements  plus  one. 

2.  y/  -  2l  veetor  eontaining  the  azimuth  angles  (in  degrees)  around  the  rotor 
disk.  This  is  a  column  vector  containing  n  elements,  where  n  equals  the  number  of 
azimuth  seetors. 

3.  V,  -  a  vector  containing  the  induced  velocity  distribution  (in  ft/sec)along 
the  rotor  blade.  This  is  a  row  veetor  eontaining  n  elements,  where  n  equals  the  number  of 
blade  elements  plus  one. 
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4.  6*  -  a  vector  containing  the  1/rev  cyclic  pitch  (in  degrees)  with  respect  to 

azimuth  angle  at  0.7  r/R.  This  is  a  column  vector  containing  n  entries,  where  n  equals  the 
number  of  azimuth  sectors. 

5-  ^  vector  containing  the  total  cyclic  pitch,  1/rev  +  2/rev  +3/rev,  (in 

degrees)  with  respect  to  azimuth  angle  at  0.7  r/R.  This  is  a  column  vector  containing  n 
entries,  where  n  equals  the  number  of  azimuth  sectors. 

6.  Pt  ~  ^  vector  containing  blade  pitch  (in  degrees)  along  the  rotor  blade. 
This  is  a  row  vector  containing  n  elements,  where  n  equals  the  number  of  blade  elements 
plus  one. 

7.  a  -  a  matrix  containing  the  angle  of  attack  (in  degrees)  distribution.  The 
matrix  is  with  m  rows  and  n  columns,  where  m  equals  the  number  of  azimuth  sectors  and 
n  equals  the  number  of  blade  elements  plus  one. 

8.  ^  vector  containing  the  total  thrust  (in  pounds)  at  each  azimuth  sector. 

This  is  a  column  vector  containing  n  elements,  where  n  equals  the  number  of  azimuth 
sectors. 

9.  -  a  vector  containing  the  total  thrust  moment  (in  ft-lbs)  at  each 

azimuth  sector.  This  is  a  column  vector  containing  n  elements,  where  n  equals  the 
number  of  azimuth  sectors. 

10.  DM^  -  a  vector  containing  the  total  drag  moment  (in  ft-lbs)  at  each 

azimuth  sector.  This  is  a  column  vector  containing  n  elements,  where  n  equals  the 
number  of  azimuth  sectors. 

11.  dT-  a  matrix  containing  the  differential  thrust  (in  pounds)  distribution. 
This  is  a  matrix  with  m  rows  and  n  columns,  where  m  equals  the  number  of  azimuth 
sectors  and  n  equals  the  number  of  blade  elements  plus  one. 

12.  dM  -  a  matrix  containing  the  differential  thrust  moment  (in  ft-lbs) 
distribution.  This  is  a  matrix  with  m  rows  and  n  columns,  where  m  equals  the  number  of 
azimuth  sectors  and  n  equals  the  number  of  blade  elements  plus  one. 
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13.  dD  -  a  matrix  containing  the  differential  drag  (in  pounds)  distribution. 
This  is  a  matrix  with  m  rows  and  n  eolumns,  where  m  equals  the  number  of  azimuth 
seetors  and  n  equals  the  number  of  blade  elements  plus  one. 

14.  Ciplot  -  a  matrix  containing  the  CL  distribution  over  the  rotor  disk.  This 
is  a  matrix  with  n  rows  and  n  eolumns,  where  n  equals  the  number  of  blade  elements. 

15.  Coplot  -  a  matrix  eontaining  the  CD  distribution  over  the  rotor  disk.  This 
is  a  matrix  with  n  rows  and  n  columns,  where  n  equals  the  number  of  blade  elements. 

When  JANRAD  2P_RVR  has  eompleted  saving  the  data,  the  Execution  Menu 
page  will  be  displayed.  See  Figure  64.  Choose  the  desired  function. 

***  EXECUTION  MENU  *** 

1.  Rotor  Performance  Analysis 

2.  Change  Data 

3.  Quit 

Enter  a  1,  2,  or  3  »| 

Figure  64.  Execution  Menu  Page. 

If  you  ehose  option  2,  you  will  be  able  to  reload  the  file  you  were  just  working  on 
or  load  an  entirely  different  file.  If  you  ehose  to  load  an  entirely  different  file,  it  must  be 
of  the  same  type  (either  Conventional  or  Compound  Helieopter)  as  the  one  you  just 
completed  working  on.  If  you  need  to  ehange  the  type  of  helieopter  (eonventional  or 
compound)  you  must  select  option  3  from  the  Exeeution  Menu  page  and  restart  JANRAD 
2P  RVR. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

The  original  intent  of  this  investigation  was  to  determine  to  what  extent  helicopter 
forward  flight  could  be  improved  using  Reverse  Velocity  Rotor  (RVR)  technology  and 
Higher  Harmonic  Stall  Control  (HHSC).  JANRAD  2P_RVR  computational  analysis 
shows  that  HHSC  alone  can  produce  increases  in  forward  flight  speeds  of  up  to  ten 
percent  in  conventional  helicopters.  Although  these  improvements  were  not  as  radical  as 
Arcidiacono  predicted  they  are  still  significant.  The  redistribution  of  the  lift  over  the 
rotor  disk  allows  the  helicopter  to  operate  with  a  greater  percentage  of  the  retreating  side 
of  the  rotor  system  engulfed  in  reverse  flow  while  still  meeting  the  rigid  trim  criteria 
imposed  by  the  program’s  original  trim  module.  Employment  of  the  2/rev  feathering 
incurred  minimal  power  penalties  and  allowed  the  blade  thrust  coefficient  to  approach 
ideal  values.  Use  of  the  HHSC  worked  equally  well  when  used  with  compound 
helicopter  configurations  up  to  175  knots,  although  the  extra  weight  and  drag  caused  by 
the  wing  increased  power  requirements. 

At  typical  helicopter  airspeeds  (less  than  160  knots)  the  RVR  is  somewhat 
inefficient.  The  RVR’s  basic  shape  does  not  perform  nearly  as  well  as  the  exotic  blades 
currently  in  use  by  rotorcraft  but,  the  use  of  a  symmetrical,  double  ended  airfoil  allows 
the  generation  of  lift  with  both  forward  and  reverse  flows  across  the  blades,  which  would 
allow  the  generation  of  lift  on  the  retreating  side  at  a  dramatically  wider  range  of 
airspeeds  than  conventionally  shaped  blades.  RVR  equipped  rotorcraft  would  require 
more  power  to  hover  and  climb  than  a  rotorcraft  equipped  with  proprietary  non-linear 
twist  or  taper.  The  true  potential  of  the  RVR  is  best  exploited  at  high  forward  flight 
speeds.  As  flight  speed  increases  the  rotor  must  be  slowed  based  on  predetermined 
optimum  advance  ratio  scheduling.  As  demonstrated  in  the  brief  JVHL  study,  the  RVR  is 
most  efficiently  placed  in  this  environment  when  it  is  used  on  an  aircraft  equipped  with  a 
compound  wing  and  auxiliary  thrust  devices.  An  RVR  equipped  aircraft  would  spend  the 
majority  of  its  time  in  cruise  flight,  minimizing  time  spent  hovering  and  at  low  airspeeds. 
At  cruise  speed  the  RVR  must  be  controlled  by  a  unique  RVR  specific  2/rev  HHSC 

control  system.  This  system  builds  upon  the  conventional  2/rev  HHSC  system  and  is 
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required  to  control  the  angle  of  attack  of  the  blades  around  the  rotor  system’s  azimuth  to 
maximize  the  RVR’s  unique  performance.  This  angle  of  attack  manipulation  allows  for 
lift  generation  in  all  quadrants  of  the  system  at  cruise  speed.  The  predictable  lift 
generation  allows  for  controlled  aerodynamic  flight  of  the  rotor  disk. 

B,  RECOMMENDATIONS 

Improvements  to  the  program  can  be  made  in  several  ways.  First,  the  graphical 
user  interface  features  within  MATLAB  should  be  harnessed  to  allow  users  to  enter  data 
via  editable  dialogue  and  text  boxes.  The  program  can  also  be  modified  to  allow  for 
offsetting  the  center  of  gravity  from  the  center  of  the  rotor  hub  during  the  trimming 
process.  As  more  data  are  made  available  more  airfoils  should  be  added  to  the  selection 
list  to  allow  for  broader  experimentation. 

Although  JANRAD  2P_RVR  has  the  capacity  to  accommodate  3/rev  inputs,  the 
study  was  limited  to  2/rev  inputs.  Investigation  into  the  effects  of  combining  3/rev  inputs 
with  both  a  non-RVR  and  RVR  system  may  yield  unforeseen  benefits  or  behaviors. 

The  development  of  a  HHSC/RVR  blade  dynamics  and  stability  and  control 
module  should  be  considered.  The  module  should  include  provisions  for  both  analytical 
and  visual  output.  These  modules  should  also  be  able  to  ‘plug-in’  to  the  existing 
JANRAD  2P_RVR  performance  module  to  give  the  user  a  complete  preliminary  design 
tool. 

With  the  advent  of  new  V/STOL  technologies  the  use  of  lift  fans,  direct  jet 
thrusters,  and  other  lift  augmentation  devices  will  become  more  commonplace.  Future 
versions  of  the  program  must  be  able  to  fully  incorporate  these  devices  and  should  be 
able  to  integrate  their  effects  with  the  traditional  helicopter  lift  (main  rotor)  and  anti¬ 
torque  (tail  rotor)  devices.  New  techniques  will  have  to  be  developed  to  combine  the 
effects  of  these  systems  to  accurately  determine  downwash  velocities,  climb 
performance,  and  controllability. 
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APPENDIX  A.  JANRAD.M 


%  Joint  Army  Navy  Rotorcraft  Analysis  and  Design 
%  JANRAD 
%  Version  2P_RVR 

%  All  modification  to  original  code  written  by 
%  MAJ  Steven  Van  Riper 
%  August  2003 

% 

%  Core  programming  technology 
%  developed  by 
%  E.  R.  Wood,  Ph.D. 

%  MAJ  Bob  Nicholoson 
%  MAJ  Walter  Wirth 
%  September  1993 
% 

%  This  program  is  an  interactive  preliminary  design  tool 
%  developed  to  aid  the  design  student  in  determination  of 
%  initial  rotorcraft  configurations  and  in  the  calculation 
%  of  performance  and  other  parameters  with  or  without 
%  IP  +  2P  +  3P  rotor  control. 

%  The  program  will  work  for  conventional  or  compound  rotorcraft. 
%  The  program  can  take  into  account  various  types  of  powered  lift 
%  including  lift  fans  and  direct  jet  thrusters. 

%  It  will  provide  accurate  data  for  airspeeds  less  than  10 
%  knots  and  greater  than  or  equal  to  50  knots. 

%  Variable  list  for  janrad.m,  trim.m,  thrcalc.m,  tmcalc.m,  dmcalc.m, 
%  hh02clcd.m,  vrl2clcd.m,  plotter,  m,  and  perf  m 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


a  lift  curve  slope  of  rotor  system  airfoil 

Adisk  area  of  rotor  disk 
Afh  fuselage  equivalent  flat  plate  drag  area 

Afv  vertical  projected  area  (fuselage  area  under  disk) 
airfoil  rotor  system  airfoil  type  (HH02WR12) 

alpha  angle  of  attack,  rotor  blade  radial  segment 
alpham  matrix  of  alpha  values  for  each  blade  element 
alphaT  rotor  tip  path  plane  angle 

alphaTfix  determines  if  user  fixes  alphaT  to  specified  value 

area  area  of  Lissajous  figures  in  plotter.m 
b  number  of  rotor  blades 

B  tip  loss  parameter 

betao  rotor  coning  angle 

betat  geometric  angle,  rotor  blade  radial  segment 

bhoriz  span,  horizontal  tail 

bvert  span,  vertical  tail 

bwing  span,  wing 

cblade  chord,  rotor  blade 

CD  drag  coefficient,  rotor  blade  radial  segment 

CDohoriz  profile  drag  coefficient,  horizontal  tail 
CDovert  profile  drag  coefficient,  vertical  tail 
CDOwing  profile  drag  coefficient,  wing 
CDhoriz  drag  coefficient,  horizontal  tail 
CDplot  matrix  of  CD  values  for  each  blade  element 
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%  CDvert  drag  coefficient,  vertical  tail 

%  CDwing  drag  coefficient,  wing 
%  CH  rotor  H  force  coefficient 

%  CH  sig  CH/solidity 

%  CL  lift  coefficient,  rotor  blade  radial  segment 

%  CLhoriz  lift  coefficient,  horizontal  tail 
%  CLplot  matrix  of  CL  values  for  each  blade  element 

%  CLvert  lift  coefficient,  vertical  tail 

%  CLwing  lift  coefficient,  wing 
%  CQ  rotor  torque  coefficient 

%  CQ_sig  CQ/solidity 

%  CT  rotor  thrust  coeeficient 

%  CT  sig  CT/solidity 

%  dD  differential  drag,  rotor  blade  radial  segment 

%  ddD  differential  drag,  rotor  blade  tip 
%  ddDM  differential  drag  moment,  rotor  blade  tip 
%  ddM  differential  thrust  moment,  rotor  blade  tip 
%  ddT  differential  thrust,  rotor  blade  tip 
%  delM  change  in  total  thrust  moment 
%  dev  AC  deviation  used  within  Adjusting  Collective  routine 

%  devCY  deviation  used  within  Adjusting  Cyclic  routine 

%  devMT  deviation  used  within  Mean  Thrust  routine 
%  devTC  deviation  used  within  Trimming  Collective  routine 
%  Dftotal  resultant  of  fuselage  drag  and  aux  thrust 
%  Dfuse  total  drag  generated  by  no  rotor  bodies 
%  DL  disk  loading 

%  dM  differentail  thrust  moment,  rotor  blade  radial  segment 

%  DMpsi  total  blade  drag  moment  at  specific  azimuth  angle 
%  dr  rotor  blade  radial  segment  width 

%  Drotor  rotor  system  drag 

%  dT  differential  thrust,  rotor  blade  radial  segment 

%  Dtplot  matrix  containing  dTs  for  all  blade  elements 

%  Dhoriz  drag,  horizontal  tail 

%  dthetadM  change  in  cyclic  pithch  with  change  in  thrust  moment 
%  Dvert  drag,  vertical  tail 

%  Dwing  drag,  wing 
%  e  effective  hinge  offset 

%  ewing  wing  efficiency  factor 

%  filename  name  of  input  file 
%  FM  figure  of  merit 

%  grip  lenght  of  inner  non-aerodynamic  protion  of  blade 
%  GW  aircraft  gross  wieght 

%  helochoice  determines  choice  of  compound  or  conventional  helo 
%  Flrotor  rotor  FI  force 

%  hub  gives  size  of  hub  in  plotter.m 

%  lamdaT  forward  flight  induced  velocity  parameter 
%  Lftotal  total  lift  generated  by  non-rotor  bodies 
%  Lhoriz  lift,  horizontal  tail 

%  lifterror  difference  between  actual  lift  and  desired  lift 
%  LoverT  Lift  over  thrust  percentage 

%  Lvert  lift,  vertical  tail 

%  Lwing  lift,  wing 

%  Mlc  first  harmonic  (cosine)  thrust  moment  coefficient 
%  Mis  first  harmonic  (sine)  thrust  moment  coefficient 
%  Machtip  Mach  number  at  rotor  blade  tip 

%  mblade  mass  of  rotor  blade 
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%  Mpsi  total  blade  thrust  moment  at  speeifie  azimuth  angle 
%  mu  advanee  ratio 

%  muehoiee  determines  if  user  would  like  JANRAD  to  ehange  mu  to 

%  a  desired  value 

%  muerror  differenee  between  aetual  mu  and  desired  mu 
%  naz  number  of  azimuth  seetors 

%  nbe  number  od  blade  elements 

%  omega  rotor  rotational  veloeity 
%  PA  pressure  altitude 

%  phi  inflow  angle,  rotor  blade  radial  segment 

%  phitip  inflow  angle,  rotor  blade  tip 

%  plotehoiee  determine  whieh  plot  to  display  in  plotter.m 

%  Protor  power  required  by  the  rotor 

%  psi  azimuth  angle 

%  q  dynamie  pressure 

%  Qrotor  rotor  torque 

%  r  radius,  rotor  blade  radial  segment 

%  R  rotor  blade  radius 

%  Rbar  Reff-e 

%  RbarT  rT*Rbar 

%  Reff  effeetive  rotor  radius,  tip  loss 

%  rho  ambient  air  density 

%  rT  loeation  of  resultant  thrust  veetor 

%  rvrCP  determines  if  normal  or  rvr  eyelie  piteh  seheduling  is  used 
%  solidity  solidity 

%  Shoriz  area,  horizontal  tail 

%  Svert  area,  vertieal  tail 

%  Swing  area,  wing 

%  T  rotor  thrust 

%  Taux  auxiliary  thrust 

%  Tauxehoiee  determines  if  user  would  like  JANRAD  to  mateh  Taux 

%  requirements 

%  temp  ambient  air  temperature 

%  theta  eyelie  piteh 

%  thetahub  hub  thetas  in  plotter.m 

%  theta  le  first  harmonie  (eosine)  of  eyelie  piteh 

%  thetals  first  harmonie  (sine)  of  eyelie  piteh 

%  theta2e  seeond  harmonie  (eosine)  of  eyelie  piteh 

%  thetaJs  third  harmonie  (sin)  of  eyelie  piteh 

%  thetaJPoffset  3/rev  offset 

%  thetao  eolleetive  piteh  at  .7  r/R 

%  thetaXP  total  vlue  of  eyelie  piteh  with  all  inputs 

%  Tpsi  total  blade  thrust  at  speeifie  azimuth  angle 

%  twist  geometrie  rotor  blade  twist 

%  Up  vertieal  eomponent  of  veloeity 

%  Uptip  vertieal  eomponent  of  veloeity  at  tip 

%  Ut  horizontal  eomponent  of  veloeity 

%  Uttip  horizontal  eomponent  of  veloeity  at  tip 

%  vertaux  vertieal  auxiliary  thrust  in  pounds 

%  vertauxehoiee  determines  if  auxiliary  vertieal  thrust  deviees  are 

%  installed  on  the  aireraft 

%  vi  indueed  veloeity 

%  Vinf  forward  airspeed 

%  Vtip  tip  speed 

%  wblade  wieght  of  rotor  blade 

%  wingehoiee  determines  if  user  wnats  to  allow  JANRAD  to  resize  the  wing 
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%  wpercent  percent  of  aircraft  gross  wieght  to  be  carried  by  the  wing 

%  clearing  all  the  variables  is  the  MATLAB  environment 

clear  all 

clc 

disp  (' ') 
disp  (' ') 
disp 


disp  (' 

* 

disp  (' 

* 

Joint  Army/Navy  Rotorcraft  Analysis  and  Design 

disp  (' 

* 

(JANRAD) 

disp  (' 

* 

disp  (' 

* 

Version  2P_RVR 

disp  (' 

* 

October  2003 

disp  (' 

* 

disp 

disp  (' ') 

% 

%  Determine  if  the  user  wants  to  analyze  conventional  or  compound 
%  helicopter  models 

% 

flag=l; 
while  flag  >  0 
disp(") 

disp('Do  you  want  to  analyze  conventional  or  compound  rotorcraft?  ') 
disp(") 

helochoice=input('  1.  conventional  or  2.  compound  »'); 

while  isempty(helochoice), 
dispC ') 

disp('You  must  enter  a  1  or  2') 

disp('Do  you  want  to  analyze  conventional  or  compound  rotorcraft?  ') 
helochoice=input('l.  conventional  or  2.  compound  »'); 
end 

if  helochoice~=l  &  helochoice~=2 
dispC ') 

dispC  ***  Enter  a  1  or  2  ***') 
else 
flag=0; 
end 
end 

clc 

checkl=l; 
while  check  1  >  0 
checkl=l; 
dispC ') 

disp('Do  you  want  to  edit  any  existing  file  or  create  a  new  one?') 

disp(") 

check=l; 

while  check  >  0 

answer=input('  1.  edit  existing  file  2.  create  new  fde  »'); 


%  ***  If  editing  an  existing  file:  get  fde  name,  display  edit 


%  menu,  allow  changes  to  selected  variables,  and  save  under 
%  desired  file  name.  Loads  to  and  saves  from  current 
%  directory  a  a  .mat  file.  *** 

if  answer==l 
clc 

dispC  ') 
disp(' ') 

dispC  ***  LOAD  INSTRUCTIONS  ***') 

dispC ') 

dispC  A.  Input  the  name  of  the  file  to  edit.') 

dispC  B.  The  file  was  saved  in  your  previous  session') 

dispC  with  a  ".mat"  extension.') 

dispC  C.  Do  not  include  the  extension  or  quaotations.') 

dispC  ') 

dispC  ex:  dsgnl') 

flag=0; 

while  flag  <  1 

filename  l=inputC  name  of  input  file:  ','s'); 
eval(['flag=existC", filename  1, '.mat");']) 
if  flag  <  1 
disp  (' ') 

disp  ('  The  file  does  not  exist,  try  again  or  <Ctrl-C>') 
disp  ('  to  exit  program.') 
end 
end 

eval(['load  ',filenamel]) 
while  check  >  0 
if  helochoice==l 
clc 

dispC ') 

dispC  *  *  *  CONVENTIONAL  HELICOPTER  EDIT  MENU  *  *  *') 

dispC ') 

dispC  1 .  pressure  altitude  2.  temperature') 

dispC  3.  airspeed  4.  gross  weight') 

dispC  5.  number  of  blades  6.  blade  radius') 

dispC  7.  blade  chord  8.  hinge  offset') 

dispC  9.  blade  grip  length  10.  blade  twist') 

dispC  11.  blade  wieght  12.  number  of  blade  elements') 

dispC  13.  rotational  velocity  14.  #  azimuth  sectors') 
dispC  15.  lift  curve  slope  16.  airfoil') 

dispC  17.  collective  pitch  18.  flatp late  area') 

dispC  19.  vert  projected  area  20.  horizontal  tail  area') 

dispC  21.  horizontal  tail  span  22.  horizontal  tail  CL') 

dispC  23 .  horizontal  tail  CDo  24.  vertical  tail  area') 

dispC  25.  vertical  tail  span  26.  vertical  tail  CL') 

dispC  27.  vertical  tail  CDo  28.  2P  input  (cosine)') 

dispC  29.  3P  input  ') 

dispC  0.  NO  CHANGES') 

alphaT  fix=0  ;itercount=0  ;helochoice= 1  ;rvrCP=2 ; 
else 
clc 

dispC ') 

dispC  ***  COMPOUND  HELICOPTER  EDIT  MENU  ***') 
dispC ') 
dispC 


1 .  pressure  altitude  2.  temperature') 
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dispC 

3.  airspeed 

4.  gross  weight') 

dispC 

5.  number  of  blades 

6.  blade  radius') 

dispC 

7.  blade  chord 

8.  hinge  offset') 

dispC 

9.  blade  grip  length 

10.  blade  twist') 

dispC 

1 1 .  blade  wieght 

12.  number  of  blade  elements') 

dispC 

13.  rotational  velocity  14.  #  azimuth  sectors') 

dispC 

15.  lift  curve  slope 

16.  airfoil') 

dispC 

17.  collective  pitch 

18.  fiatplate  area') 

dispC 

19.  vert  projected  area  20.  horizontal  tail  area') 

dispC 

21.  horizontal  tail  span  22.  horizontal  tail  CL') 

dispC 

23.  horizontal  tail  CDo  24.  vertical  tail  area') 

dispC 

25.  vertical  tail  span 

26.  vertical  tail  CL') 

dispC 

27.  vertical  tail  CDo 

28.  2P  input') 

dispC 

29.  3P  input 

30.  auxiliary  thrust') 

dispC 

3 1 .  wing  area 

32.  wing  span') 

dispC 

33.  wing  CL 

34.  wing  CDo') 

dispC 

35.  wing  efficiency  factor') 

dispC 

0.  NO  CHANGES') 

alphaT  fix=0  ;itercouiit=0  ;helochoice=2  ;rvrCP=2 ; 
end 

choice=input('  Select  the  parameter  to  change: '); 
if  isempty(choice), 
choice=0; 
end 

ifchoice==l, 

clc 

dispC ') 

PA 

tmp=PA; 

PA=input('Pressure  altitude  (ft); '); 
if  isempty(PA), 

PA=tmp; 

end 

elseif  choice==2, 
clc 

dispC ') 
temp 

tmp=temp; 

temp=input('temperature  (deg  F):'); 
if  isempty(temp)„ 
temp=tmp; 
end 

elseif  choice==3, 
clc 

dispC ') 

Vinf=Vinf/1.69 

tmp=Vinf; 

Vinf=input('Airspeed  (knots):  ')*1.69; 
if  isempty(Vinf), 

Vinf=tmp*1.69; 

end 

elseif  choice==4, 
clc 

dispC ') 

GW 

tmp=GW; 
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GW=input('Aircraft  Groos  weight  (lbs):  '); 
if  isempty  (GW), 

GW=tmp; 

end 

elseif  choice==5, 
cle 

dispC ') 
b 

tmp=b; 

b=input('Number  of  blade: '); 
if  isempty(b), 
b=tmp; 
end 

elseif  choice==6, 
cle 

dispC ') 

R 

tmp=R; 

R=input('Blade  radius;  center  of  hub  to  blade  tip  (ft):  '); 
if  isempty(R), 

R=tmp; 

end 

elseif  choice==7, 
cle 

dispC ') 

cblade 

tmp=cblade; 

cblade=input('Blade  chord  (ft); '); 
if  isempty(cblade), 
cblade=tmp; 
end 

elseif  choice==8, 
cle 

dispC ') 
e 

tmp=e; 

e=inputCHinge  offset  (ft); '); 
if  isempty(e), 
e=tmp; 
end 

elseif  choice==9, 
cle 

dispC ') 
grip 

tmp=grip; 

grip=inputCnon-aerodyn  inboard  portion  of  blade  (ft):  '); 
if  isempty  (grip), 
grip=tmp; 
end 

if  grip  <  le-1, 
grip=le-10; 
end 

elseif  choice==10 
cle 

dispC ') 

twist=-twist*57.3 
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tmp=twist; 

twist=input('blade  twist  (deg): '); 
if  isempty  (twist), 
twist=abs(tmp)/ 57.3; 
else 

twist=abs(twist)/57.3 ; 
end 

elseif  ehoiee==ll, 
ele 

dispC ') 
wblade 
tmp=wblade; 

wblade=input('Weight  of  aero  portion  of  blade  (lbs): '); 
if  isempty(wblade), 
wblade=tmp; 
end 

elseif  ehoice==12, 
ele 

dispC ') 
nbe 

tmp=nbe; 

nbe=input('number  of  blade  elements: '); 
if  isempty(nbe), 
nbe=tmp; 
end 

elseif  ehoiee==13, 
ele 

dispC ') 
omega 
tmp=omega; 

omega=mput('Rotor  rotational  veloeity  (rad/see):  '); 
if  isempty  (omega), 
omega=tmp; 
end 

elseif  elioiee==14, 
ele 

dispC ') 
naz 

tmp=naz; 

naz=inputCnumber  of  azimuth  seetors:  '); 
if  isempty  (naz), 
naz=tmp; 
end 

elseif  ehoice==l  5, 
ele 

dispC ') 
a 

tmp=a; 

a=inputCAverage  lift  eurve  slope  (CL  vs  alpha); '); 
if  isempty(a), 
a=tmp; 
end 

elseif  ehoiee==l  6, 
ele 

disp(") 

airfoil 
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tmp=airfoil; 
flag=l; 
while  flag  >  0 

airfoil=input('Airfoil  1.  HH-02  2.  VR-12  3.  SC1095  4.  0012  »'); 
if  isempty(airfoil), 
airfoil=tmp; 
end 

if  airfoil==l, 
flag=0; 

elseif  airfoil==2, 
flag=0; 

elseif  airfoil==3, 
flag=0; 

elseif  airfoil==4, 
flag=0; 
else 
dispC ') 

dispC  ***  Enter  a  1,  2,  3,  or  4  ***  ') 
end 
end 

elseif  choice==16  &  helochoice==2, 
cle 

disp(") 
airfoil 
tmp=airfoil; 
flag=l; 
while  flag  >  0 

airfoil=input('Airfoil  1.  HH-02  2.  VR-12  3.  SC1095  4.  0012  5.  RVR»'); 
if  isempty(airfoil), 
airfoil=tmp; 
end 

if  airfoil==l, 
flag=0; 

elseif  airfoil==2, 
flag=0; 

elseif  airfoil==3, 
flag=0; 

elseif  airfoil==4, 
flag=0; 

elseif  airfoil==5, 
flag=0; 
else 
dispC ') 

dispC  ***  Enter  a  1,  2,  3,  4,  or  5  ***  ') 
end 
end 

elseif  choice==l  7, 
cle 

dispC ') 

thetao=thetao  *57.3 
tmp=thetao; 

thetao=inputCeollective  pitch  at  .7  r/R  (deg):  ')/57.3; 
if  isempty(thetao), 
thetao=tmp/57.3; 
end 

elseif  choice==18, 
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clc 

dispC ') 

Afh 

tmp=Afh; 

Afh=input('Aircraft  equivalent  flatplate  area  (ft^2):'); 
if  isempty(Afh), 

Afli=tmp; 

end 

elseif  ehoiee==19, 
ele 

dispC ') 

Afv 

tmp=Afv; 

Afv=input('Verieal  projeeted  area  (ft^2):'); 
if  isempty(Afv), 

Afv=tmp; 

end 

elseif  ehoiee==20, 
ele 

dispC ') 

Shoriz 

tmp^Shoriz; 

Shoriz=input('Horizontal  tail  area  (ft^2):  '); 
if  isempty(Shoriz), 

Shoriz=tmp; 

end 

if  Shoriz  <  le-10, 

Shoriz=le-10; 

end 

elseif  ehoiee==21, 
ele 

dispC ') 

bhoriz 

tmp=bhoriz; 

bhoriz=inputCHorizontal  tail  span  (ft):  '); 
if  isempty(bhoriz), 
bhoriz=tmp; 
end 

if  bhoriz  <  le-10, 
bhoriz=  le-10; 
end 

elseif  ehoiee==22, 
ele 

dispC ') 

CLhoriz 

tmp=CLhoriz; 

CLhoriz=inputCExpeeted  CL  for  the  horizontal  tail:  '); 
if  isempty  (CLhoriz), 

CLhoriz=tmp; 

end 

elseif  ehoiee==23, 
ele 

dispC ') 

CDohoriz 

tmp=CDohoriz; 

CDohoriz=input('Horizontal  tail  prolfde  drag  eoef  (CDo):  '); 
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if  isempty(CDohoriz), 

CDOhoriz=tmp; 

end 

elseif  choice==24, 
clc 

dispC ') 

Svert 

tmp=Svert; 

Svert=input('Vertical  tail  area  (ft^2):  '); 
if  isempty(Svert), 

Svert=tmp; 

end 

if  Svert  <  le-10, 

Svert=le-10; 

end 

elseif  choice==25, 
clc 

dispC ') 
bvert 

tmp=bvert; 

bvert=input(' Vertical  tail  span  (ft):  '); 
if  isempty(bvert), 
bvert=tmp; 
end 

if  bvert  <  le-10, 
bvert=  le-10; 
end 

elseif  choice==26 
clc 

dispC ') 

CLvert 

tmp=CLvert; 

CLvert=input('Expected  CL  for  the  vertical  tail:  '); 
if  isempty(CLvert), 

CLvert=tmp; 

end 

elseif  choice==27, 
clc 

dispC ') 

CDovert 

tmp=CDovert; 

CDovert=inputC Vertical  tail  profile  drag  coef  (CDo); '); 
if  isempty(CDovert), 

CDovert=tmp; 

end 

elseif  choice==28, 
clc 

dispC ') 

twoPinput 

tmp=twoPinput; 

twoPinput=inputCLongitudinal  2P  harmonic  input  (degs):  '); 
if  isempty  (twoPinput), 
twoPinput=tmp; 
end 

elseif  choice==29 
clc 
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disp(") 

threePinput 

tmp=threePinput; 

threePinput=input('3P  harmonic  input  (degs):  '); 

disp(") 

threePoffset 

tmp=threePoffset; 

threePoffset=input('3P  input  offset,  must  be  between  0  and  119  degs: '); 
if  isempty  (threePinput) 
threePinput=tmp ; 
end 

elseif  helochoice==2  &  choice==30, 
clc 

dispC ') 

Taux 

tmp=Taux; 

Taux=input('Auxiliary  thrust  (lbs):  '); 
if  isempty  (Taux), 

Taux=tmp; 

end 

elseif  helochoice==2  &  choice==31, 
clc 

dispC ') 

Swing 

tmp=Swing; 

Swing=input('Wing  Area  (ft^2):'); 
if  isempty  (Swing), 

Swing=tmp; 

end 

if  Swing  <  le-10 
Swing=le-10; 
end 

elseif  helochoice==2  &  choice==32, 
clc 

dispC ') 

bwing 

tmp=bwing; 

bwing=input('Wing  span  (ft); '); 
if  isempty(bwing), 
bwing=tmp; 
end 

ifbwing  <  le-10, 
bwing=le-10; 
end 

elseif  helochoice==2  &  choice==33, 
clc 

dispC ') 

CLwing 

tmp=CLwing; 

CLwing=inputCexpected  CL  for  the  wing:'); 
if  isempty(CLwing), 

CLwing=tmp; 

end 

elseif  helochoice==2  &  choice==34, 
clc 

dispC ') 


74 


CDowing 

tmp=CDowing; 

CDowing=input('Wing  profile  drag  coef  (CDo):'); 
if  isempty(CDowing), 

CDowing=tmp; 

end 

elseif  helochoice==2  &  choice==35, 
clc 

dispC ') 

ewing 

tmp=ewing; 

ewing=input('Wing  efficiency  factor  (e): '); 
if  isempty(ewing), 
ewing=tmp; 
end 

if  ewing  <  le-10, 
ewing=le-10; 
end 

elseif  helochoice==2  &  choice==0 
GW1=GW; 
clc 

ifVinf>  16.9 
disp(") 

disp('Do  you  want  to  fix  the  Tip  Path  Plane  (recommended  for  compound 

helicopters),') 

alphaTfix=rnput('l.  yes  or  2.  no?'); 
if  alphaTfix==l 
disp(") 

alphaTinput=rnput('Enter  the  Tip  Path  Plane  angle  in  degrees:  '); 
alphaTinput=alphaTinput/57.3 ; 

dispC***  You  have  fixed  the  Tip  Path  Plane  angle,  please  set  auxiliary  ***') 
dispC***  thrust  to  equal  total  drag  to  avoid  erroneous  solutions.  ***1^ 
else 
disp(") 

disp('Tip  path  plane  will  be  determined  by  JANRAD') 
disp('Press  any  key  to  continue...') 
pause 
end 

%  added  14-16  Sep  2003  to  allow  user  to  disable  the  retrim  routine,  select  Taux  based 
%  on  drag,  and  allow  for  RVR  based  cyclic  inputs  in  trim.m 
flag=l; 
while  flag  >  0 
disp(") 

disp('Do  you  want  JANRAD  to  set  auxiliary  thmst  to  equal  total  drag, ') 
Tauxchoice=input('l.  yes  or  2.  no?'); 
while  isempty(Tauxchoice), 
dispC ') 

disp('You  must  enter  a  1  or  2') 

disp('Do  you  want  JANRAD  to  set  auxiliary  thrust  to  equal  total  drag, ') 
Tauxchoice=input('l.  yes  or  2.  no?'); 
end 

if  Tauxchoice~=l  &  Tauxchoice~=2 
dispC ') 

dispC  ***  Enter  a  1  or  2  ***') 
elseif  T  auxchoice==  1 

dispCJANRAD  will  adjust  auxiliary  thrust  to  equal  total  drag') 
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disp('Press  any  key  to  continue...') 

pause 

flag^O; 

elseif  Tauxchoice==2 

dispC***  Auxiliary  thrust  defined  by  user.  ***') 
disp('Press  any  key  to  eontinue...') 
pause 
flag=0; 
else 
flag^O; 
end 
end 

rvrCP=2; 
ifVinf/1.69>200 
rvrCP^O; 
flag=l; 
while  flag  >  0 
disp(") 

disp('Do  you  want  JANRAD  to  use  RVR  based  (eonstrained)  2/rev  Cyelie 

Piteh, ') 

rvrCP=input('l.  yes  or  2.  no?'); 
while  isempty(rvrCP), 
dispC  ') 

disp('You  must  enter  a  1  or  2') 

disp('Do  you  want  JANRAD  to  use  RVR  based  (eonstrained)  2/rev  Cyelie 

Piteh, ') 

rvrCP=input('l.  yes  or  2.  no?'); 
end 

if  rvrCP~=l  &  rvrCP~=2 
dispC ') 

dispC  ***  Enter  a  1  or  2  ***') 
elseifrvrCP==l 
disp(") 

dispCRVR  2/rev  Cyelie  Piteh  will  be  used') 
dispCPress  any  key  to  eontinue...') 
pause 
flag=0; 

elseif  rvrCP==2 
dispC) 

dispCConventional  2/rev  eyelie  piteh  will  be  used') 
dispCPress  any  key  to  eontinue...') 
pause 
flag=0; 
else 
flag=0; 
end 
end 
else 
end 
else 
end 
ele 

dispC ') 
dispC ') 

dispC  ***  SAVE  INSTRUCTIONS  ***') 

dispC ') 
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dispC  A.  Save  the  new  data  to  a  specificed  file  name.') 

dispC  B.  Do  not  use  an  extension  or  quotations.') 

dispC  C.  Use  the  letter/number  combinations  of  20  characters  or  less.') 

dispC  D.  The  file  will  be  saved  with  a  ".mat"  extension.') 

dispC  ') 

dispC  ex:  dsgn_2') 
dispC ') 

dispC  E.  If  you  made  no  changes,  press  <  Enter  >,  the  file  will') 
dispC  be  saved  with  the  original  name.') 

flag=l; 
while  flag  >  0 

filenameO=filename  1 ; 
filename l^inputC  save  file  as:  ','s'); 
if  isempty(filenamel) 
filename  1  =filenameO ; 
end 

clear  filenameO 
if  length(filenamel)  >  20, 
dispC  ') 

dispC  use  20  characters  or  less') 
flag=l; 
else 
flag=0; 
end 
end 

eval(['save  ',filenamel]) 
check=0; 
elseif  choice==0 
if  helochoice==l 
Swing=le-10; 
bwing=le-10; 

CLwing=0; 

CDowing=0; 

ewing=le-10; 

Taux=0; 

GW1=GW; 

else 

end 

clc 

dispC ') 
dispC ') 

dispC  ***  SAVE  INSTRUCTIONS  ***') 

dispC ') 

dispC  A.  Save  the  new  data  to  a  specificed  file  name.') 

dispC  B.  Do  not  use  an  extension  or  quotations.') 

dispC  C.  Use  the  letter/number  combinations  of  20  characters  or  less.') 

dispC  D.  The  file  will  be  saved  with  a  ".mat"  extension.') 

dispC  ') 

dispC  ex:  dsgn_2') 
dispC ') 

dispC  E.  If  you  made  no  changes,  press  <  Enter  >,  the  file  will') 
dispC  be  saved  with  the  original  name.') 

flag=l; 
while  flag  >  0 

filenameO=filename  1 ; 

filename l=inputC  save  file  as:  ','s'); 
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if  isempty(filenamel) 
filename  1  =filenameO ; 
end 

clear  filenameO 
if  length(filenamel)  >  20, 
dispC  ') 

dispC  use  20  characters  or  less') 
flag=l; 
else 
flag=0; 
end 
end 

eval(['save  ',filenamel]) 
check=0; 

else 
dispC ') 

dispC  Enter  a  displayed  number  ...press  any  key  to  continue') 
pause 
end 
end 


%  ***  If  creating  a  new  file:  get  input  for  required  variables 
%  and  save  under  desired  file  name.  Save  to  current 
%  directory  as  a  .mat  file.  *** 

elseif  answer==2, 

alphaT  fix=0  ;itercount=0 ; 
clc 

PA=inputCPressure  Altitude  (ft):  '); 
while  isempty(PA), 
dispC  ') 

dispCYou  must  enter  a  numerical  value') 
PA=inputCPressure  altitude  (ft):  '); 
end 

temp=input('Temperature  (deg  F):'); 
while  isempty(temp), 
dispC  ') 

dispCYou  must  enter  a  numerical  value') 
temp=input('Temperature  (deg  F):  '); 
end 

Vinf=input('Airspeed  (knots):  ')*1.69; 
while  isempty(Vinf), 
dispC  ') 

dispCYou  must  enter  a  numerical  value') 
Vinf=input('Airspeed  (knots):  ')*1.69; 
end 

GW=input('Aircraft  gross  wieght  (lbs):  '); 
while  isempty(GW), 
dispC  ') 

dispCYou  must  enter  a  numerical  value') 
GW=input('Aircrafl  gross  wieght  (lbs):'); 
end 

b=input('Number  of  Blades:  '); 
while  isempty(b), 
dispC ') 
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disp('You  must  enter  a  numerical  value') 
b=input('Number  of  Blades: '); 
end 

R=input('Blade  radius:  center  of  hub  to  blade  tip  (ft):  '); 
while  isempty(R), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 

R=input('Blade  radius:  center  of  hub  to  blade  tip  (ft):  '); 
end 

cblade=input('Blade  chord  (ft):  '); 
while  isempty(cblade), 
dispC ') 

disp('You  must  enter  a  numerical  value') 
cblade=input('Blade  chord  (ft):  '); 
end 

e=input('Hinge  offset  (ft)'); 
while  isempty(e) 
dispC ') 

disp('You  must  enter  a  numercial  value') 
e=input('Hinge  offset  (ft)'); 
end 

grip=input('Non-aerodynamic  inboard  portion  of  blade  (ft):  '); 
while  isempty(grip) 
dispC  ') 

dispCYou  must  enter  a  numercial  value') 
grip=input('Non-aerodynamic  inboard  portion  of  blade  (ft):  '); 
end 

while  grip  <  le-10, 
grip=le-10; 
end 

twist=input('Blade  twist  (deg):  ')/57.3;,twist=abs(twist); 
while  isempty(twist), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 
twist=input('Blade  twist  (deg):  ')/57.3;,twist=abs(twist); 
end 

wblade=input('weight  of  aero  portion  of  one  blade  (lbs): '); 
while  isempty(wblade), 
dispC ') 

dispCYou  must  enter  a  numerical  value') 
wblade=input('weight  of  aero  portion  if  one  blade  (lbs): '); 
end 

nbe=input('Number  of  blade  elements:  '); 
while  isempty(nbe), 
dispC ') 

dispCYou  must  enter  a  numerical  value') 
nbe=input('Number  of  blade  elements: '); 
end 

omega=input('Rotor  rotational  velocity  (rad/sec): '); 
while  isempty(omega), 
dispC  ') 

dispCYou  must  enter  a  numerical  value') 
omega=input('Rotor  rotational  velocity  (rad/sec):  '); 
end 

naz=input('Number  of  azimuth  sectors:  '); 
while  isempty(naz), 
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dispC ') 

disp('You  must  enter  a  numerical  value') 
naz=input('Number  of  azimuth  sectors:  '); 
end 

a=input('Lift  curve  slope  of  rotor  airfoil  (CL  vs  alpha):  '); 
while  isempty(a), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 
a=input('Lift  curve  slope  of  rotor  airfoil  (CL  vs  alpha):  '); 
end 

flag=l; 
while  flag  >  0 
if  helochoice==l 

airfoil=input('Airfoil  1.  HH-02  2.  VR-12  3.  SC1095  4.  0012  »'); 
while  isempty(airfoil), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 
airfoil=input('Airfoil  1.  HH-02  2.  VR-12  3.  SC1095  4.  0012  »'); 
end 

ifairfoil~=l  &  airfoil~=2  &  airfoil~=3  &  airfoil~=4 
dispC  ') 

dispC  ***  Enter  a  1,  2,  3,  or  4  ***') 
else 
flag=0; 
end 
else 

airfoil=input('Airfoil  1.  HH-02  2.  VR-12  3.  SC1095  4.  0012  5.  RVR  »'); 
while  isempty(airfoil), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 

airfoil=input('Airfoil  1.  HH-02  2.  VR-12  3.  SC1095  4.  0012  5.  RVR  »'); 
end 

if  airfoil~=l  &  airfoil~=2  &  airfoil~=3  &  airfoil~=4  &  airfoil~=5 
dispC  ') 

dispC  ***  Enter  a  1,  2,  3,  4,  or  5  ***') 
else 
flag=0; 
end 
end 
end 

thetao=input('collective  pitch  at  .7  r/R  (deg):  ')/57.3; 
while  isempty(thetao), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 
thetao=input('Collective  pitch  at  .7  r/R  (deg):  ')/57.3; 
end 

Afh=input('Aircraft  equivalent  flatplate  area  (fC2):'); 
while  isempty(Afh) 
dispC  ') 

disp('You  must  enter  a  numerical  value') 

Afh=input('Aircraft  equivalent  flatplate  area  (ft^2):'); 
end 

Afv=input(' Vertical  projected  area  (fC2):  '); 
while  isempty(Afv) 
dispC  ') 

disp('Tou  must  enter  a  numerical  value') 
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Afv=input(' Vertical  projected  area  (ft^2):  '); 
end 

Shoriz=input('Horizontal  tail  area,  0  if  none  (ft^2): '); 
while  isempty(Shoriz), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 
Shoriz=input('Horizontal  tail  area,  0  if  none  (ft^2):  '); 
end 

if  Shoriz~=0, 

bhoriz=input('Horizontal  tail  span  (ft):  '); 
while  isempty(bhoriz), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 
bhoriz=input('Horizontal  tail  span  (ft):  '); 
end 

if  bhoriz  <  le-10, 
bhoriz=le-10; 
end 

CLhoriz=input('Expected  CL  for  the  horizontal  tail:  '); 
while  isempty(CLhoriz), 
dispC  ') 

dispC  You  must  enter  a  numerical  value') 
CLhoriz=input('Expected  CL  for  the  horizontal  tail:  '); 
end 

CDohoriz=input('Horizontal  tail  profile  drag  coef  (CDo):  '); 
while  isempty(CDohoriz), 
dispC  ') 

dispC  You  must  enter  a  numerical  value') 
CDohoriz=input('Horizontal  tail  profde  drag  coef  (CDo):  '); 
end 
else 

Shoriz=  le-10; 
bhoriz=le-10; 

CLhoriz=0; 

CDohoriz=0; 

end 

Svert=input(' Vertical  tail  area,  0  if  none  (fC2):  '); 
while  isempty(Svert), 
dispC  ') 

dispCYou  must  enter  a  numerical  value') 

Svert=input(' Vertical  tail  area,  0  if  none  (ft^2):  '); 
end 

if  Svert~=0 

bvert=input(' Vertical  tail  span  (ft):  '); 
while  isempty(bvert), 
dispC  ') 

dispC  You  must  enter  a  numerical  value') 
bvert=input(' Vertical  tail  span  (ft):  '); 
end 

if  bvert  <  le-10, 
bvert=le-10; 
end 

CLvert=input('Expected  CL  for  the  vertical  tail:  '); 
while  isempty(CLvert), 
dispC  ') 

dispC  You  must  enter  a  numerical  value') 
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CLvert=input('Expected  CL  for  the  vertical  tail:  '); 
end 

CDovert=input(' Vertical  tail  profde  drag  coef  (CDo):'); 
while  isempty(CDovert), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 

CDovert=input(' Vertical  tail  profde  drag  coef  (CDo):'); 
end 
else 

Svert=le-10; 

bvert=le-10; 

CLvert=0; 

CDovert=0; 

end 

twoPinput=input('2P  longitudinal  input  (degs):  '); 
while  isempty(twoPinput), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 
twoPinput=input('2P  longitudinal  input  (degs):  '); 
end 

threePinput=input('3P  input  (degs): '); 
while  isempty(threePinput), 
dispC  ') 

disp('You  must  enter  a  numerical  value') 
threePinput=input('3P  input  (degs): '); 
end 

threePoffset=input('3P  offset,  must  be  between  0  and  119  degs: '); 
while  isempty(threePoffset)  or  threePoffset>119  or  threePoffseKO, 
disp(") 

dispC  You  must  enter  a  numerical  value  between  0  and  119  degs') 
threePoffset=input('3P  offset  (degs): '); 
end 

Swing=le-10; 

bwing=le-10; 

CLwing=0; 

CDowing=0; 

ewing=le-10; 

Taux=0; 

rvrCP=2; 

if  helochoice==2 

Taux=input('Auxiliary  thrust  (lbs):  '); 
while  isempty(Taux), 
dispC  ') 

dispC  You  must  enter  a  numerical  value') 
Taux=input('Auxiliary  thrust  (lbs):  '); 
end 

Swing=input('Wing  area,  0  if  no  wing  (ft^2):  '); 
while  isempty(Swing) 
dispC  ') 

dispC  You  must  enter  a  numerical  value') 

Swing=input('Wing  area,  0  if  no  wing  (fC2):  '); 
end 

if  Swing~=0, 

bwing=input('Wing  span  (ft):'); 
while  isempty(bwing), 
dispC  ') 
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disp('You  must  enter  a  numerieal  value') 
bwing=input('Wing  span  (ft): '); 
end 

ifbwing  <  le-10, 
bwing=le-10; 
end 

CLwing=input('Expeeted  CL  for  the  wing:  '); 
while  isempty(CLwing), 
dispC  ') 

disp('You  must  enter  a  numerieal  value') 

CLwing=input('Expeeted  CL  for  the  wing:  '); 
end 

CDowing=input('Wing  profde  drag  eoef  (CDo):  '); 
while  isempty(CDowing), 
dispC  ') 

disp('You  must  enter  a  numerieal  value') 

CDo wing=input(' Wing  profile  drag  eoef  (CDo):  '); 
end 

ewing=input('Wing  effieieney  faetor  (e): '); 
while  isempty(ewing), 
dispC  ') 

disp('You  must  enter  a  numerieal  value') 
ewing=input('Wing  effieieney  faetor  (e): '); 
end 

if  ewing  <  le-10, 
ewing=le-10; 
end 
else 

Swing=le-10; 

bwing=le-10; 

CLwing=0; 

CDowing=0; 
ewing=  le-10; 
end 

GW1=GW; 

ele 

ifVinf>  16.9 
disp(") 

disp('Do  you  want  to  fix  the  Tip  Path  Plane  (recommended  for  compound 

helicopters),') 

alphaTfix=input('l.  yes  or  2.  no?'); 
if  alphaTfix==l 
disp(") 

alphaTinput=mput('Enter  the  Tip  Path  Plane  angle  in  degrees:  '); 
alphaTinput=alphaTinput/5  7 . 3 ; 

dispC***  You  have  fixed  the  Tip  Path  Plane  angle,  please  set  auxiliary  ***') 
dispC***  thrust  to  equal  total  drag  to  avoid  erroneous  solutions.  ***1^ 
else 
disp(") 

disp('Tip  path  plane  will  be  determined  by  JANRAD') 
disp('Press  any  key  to  continue...') 
pause 
end 

%  added  14-16  Sep  2003  to  allow  user  to  disable  the  retrim  routine,  select  Taux  based 
%  on  drag,  and  allow  for  RVR  based  cyclic  inputs  in  trim.m 
flag=l; 
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while  flag  >  0 
disp(") 

disp('Do  you  want  JANRAD  to  set  auxiliary  thrust  to  equal  total  drag, ') 
Tauxchoiee=input('l.  yes  or  2.  no?'); 
while  isempty(Tauxchoiee), 
dispC  ') 

disp('You  must  enter  a  1  or  2') 

disp('Do  you  want  JANRAD  to  set  auxiliary  thrust  to  equal  total  drag, ') 
Tauxchoice=input('l.  yes  or  2.  no?'); 
end 

if  Tauxchoice~=l  &  Tauxchoice~=2 
dispC  ') 

dispC  ***  Enter  a  1  or  2  ***') 
elseif  Tauxchoice==l 

disp('JANRAD  will  adjust  auxiliary  thmst  to  equal  total  drag') 

disp('Press  any  key  to  continue...') 

pause 

flag=0; 

elseif  Tauxchoice==2 

dispC***  Auxiliary  thrust  defined  by  user.  ***') 
disp('Press  any  key  to  continue...') 
pause 
flag=0; 
else 
flag=0; 
end 
end 

rvrCP=2; 
ifVinP1.69>200 
rvrCP^O; 
flag=l; 
while  flag  >  0 
disp(") 

disp('Do  you  want  JANRAD  to  use  RVR  based  (constrained)  2/rev  Cyclic  Pitch, 

') 

rvrCP=input('l.  yes  or  2.  no?'); 
while  isempty(rvrCP), 
dispC ') 

dispCYou  must  enter  a  1  or  2') 

dispCDo  you  want  JANRAD  to  use  RVR  based  (constrained)  2/rev  Cyclic 

Pitch, ') 

rvrCP=input('l.  yes  or  2.  no?'); 
end 

if  rvrCP~=l  &  rvrCP~=2 
dispC ') 

dispC  ***  Enter  a  1  or  2  ***') 
elseif  rvrCP==l 
disp(") 

disp('RVR  2/rev  Cyclic  Pitch  will  be  used') 
disp('Press  any  key  to  continue...') 
pause 
flag=0; 

elseif  rvrCP==2 
disp(") 

disp('Conventional  2/rev  cyclic  pitch  will  be  used') 
disp('Press  any  key  to  continue...') 
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pause 

flag^O; 

else 

flag^O; 

end 

end 

else 

end 

else 

end 

else 

end 

ele 

dispC ') 
dispC ') 

dispC  ***  Save  Instruetions  ***') 

dispC ') 

dispC  A.  Save  the  data  to  a  speeified  file  name.') 

dispC  B.  Do  not  use  an  extension  or  qoutations.') 

dispC  C.  Use  letter/number  eombinations  of  6  eharaeters  or  less.') 

dispC  D.  The  file  will  be  saved  with  a  ".mat"  extension.') 

dispC ') 

dispC  ex:  dsgn  T) 
dispC ') 

dispC  E.  If  you  do  not  enter  a  name,  the  default  is  "dsgn  l '") 

flag=l; 

while  flag  >  0 

filename l=inputC  save  file  as:  ','s'); 
if  isempty(filenamel) 
filename  1  ='dsgn_  1 '; 
end 

if  length(filenamel)  >  20, 
disp(") 

dispCUse  20  eharaeters  or  less.') 
flag=l; 
else 
flag=0; 
end 
end 

eval(['save  ',filenamel]) 
eheek=0; 
else 
dispC ') 

dispC  Enter  a  1  or  2') 
end 
end 
ele 

dispC ') 
pause(3) 
eheek2=l; 
while  eheek2  >  0 
ele 

dispC ') 

dispC  ***  EXECUTION  MENU  ***') 

dispC ') 

dispC  1 .  Rotor  Performanee  Analysis') 
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dispC  2.  Change  Data') 

disp('  3.  Quit') 

check3=l; 
while  eheek3  >  0 

answer=input('  Enter  a  1,  2,  or  3  »'); 
save  temp  eheekl  eheek2  eheek3  fdenamel 
if  answer==l, 
perf 

load  temp 
eheek3=0; 
elseif  answer==2 
ele 

eheek2=0; 
eheek3=0; 
elseif  answer==3, 
dispC ') 

dispC  Thank  you  for  using  JANRAD') 
eheekl=0; 
eheek2=0; 
eheek3=0; 
end 
end 
end 
end 
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APPENDIX  B.  DMCALC.M 


%  dmcalc  calculates  the  total  drag  along  a  blade  at 
%  each  azimuth  (psi)  location 

Up=zeros(size(psi  *r)) ; 

Ut=zeros(size(Up)); 

alpham=zeros(size(Up)); 

%  added  for  plotting  routine 
alphamplot=zeros(size(Ut)); 

dD=zeros(size(Up)); 

ddD=zeros(size(psi)); 

ddDM=zeros(size(psi)); 

for  i=l:length(psi), 

Up(i,:)=vi.*cos(betao)+Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Ut(i,:)=r.*omega+Vinf*cos(alphaT)*sin(psi(i)); 

mach  =  ( Vtip .  *cos(alphaT) .  *r./R+ V inf  *  sin(psi(i)))/(49 . 05  *  sqrt(temp+460)) ; 
phi=atan2(Up(i, :  ),Ut(i, : )) ; 

%  Replace  IP  theta  with  2P  theta 
%  alpha=theta(i)+betat-phi; 
alpha=thetaXP(i)+b  etat-phi ; 

%  Test  for  plotting  routine 
alpha  1  =thetaXP(i)+betat-phi; 
alphamplot(i,  :)=alphal ; 

alpham(i,  :)=alpha; 


ifairfoil==l, 

[CL,CD]=hh02clcd(alpha); 

%  additional  airfoils  added  for  experimentation 
elseif  airfoil==2, 

[CL,CD]=vrl2clcd(alpha); 
elseif  airfoil==3, 

[CL,CD]=scl095r8clcd(alpha,mach); 
elseif  airfoil==4, 

[CL,CD]=oo  1 2clcd(alpha,mach); 
elseif  airfoil==5, 

[CL,CD]=rvrclcd(alpha,mach); 

end 

%  Populate  CLplot  and  CDplot  matrix 
CLplot(i,:)=CL; 

CDplot(i,:)=CD; 

dD(i,:)=0.5*rho*cblade*dr*(Up(i,:)  *2).*(CL.*sin(phi)+CD.*cos(phi)); 
dDM=dD(i,:).*r; 

DMpsi(i)=sum(dDM) ; 

%  ***  calculations  for  tip  loss  area  *** 
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Uptip=Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Uttip=(R-(r-Reff)/2)*omega+Vinf*cos(alpbaT)*sin(psi(i)); 

pbitip=atan2(Uptip,Uttip); 

ddD(i)=0.5*rbo*cblade*(R-Reff)*(Uptip''2+Uttip.^2)*(.009*cos(pbitip')); 

ddDM(i)=ddD(i)*(R-(R-Reff)/2); 

DMpsi(i)=DMpsi(i)+ddDM(i); 

end 
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APPENDIX  C.  PERF.M 


%  ***  Performance  output  subroutine  *** 

%  ***  Clearing  all  variable  in  the  MATLAB  environment  *** 
clear 

%  ***  Loading  input  and  check  parameters  *** 
load  temp 

eval(['load  ',fdenamel]) 
trim 

disp('trimming  complete') 

%  ***  Calculation  of  output  parameters  *** 

Protor=mean(D]Vlpsi)  *b  *omega/ 550; 

Qrotor=mean(DMpsi)  *b ; 
sohdity=b*cblade/(pi*R); 

CQ=Qrotor/(Adisk*rho  *  V  tip*2  *R) ; 

CH=Hrotor/(Adisk*rho*Vtip'^2); 

CT_sig=CT/  solidity; 

CQ_sig=CQ/ solidity; 

CH_sig=CH/ solidity; 

Machtip=(Vtip  *cos(alphaT)+V  inf)/(49 .05  *  sqrt(temp+460)) ; 
ifVinf<  16.9, 

DL=T/(pi*R^2); 

FM=(T  *  sqrt(DL/(2  *rho)))/(5 50 *Protor); 
else 
DL=0; 

FM=0; 

end 

clc 

dispC ') 

dispC  ***  PERFORMANCE  CALCULATIONS  COMPLETE  ***') 
disp(' ') 

dispC  Do  you  want  the  performance  results  displayed  on  screen?') 

flag=l; 

while  flag  >  0 

answer=input('  1.  yes  2.  no  »'); 
if  answer==2, 
flag=0; 

elseif  answer==l 
%  ***  output  to  screen  *** 


clc 

dispC ') 

dispC  ***  INPUT  DATA  ***') 
eval(['disp("  ',fdename  1 ,'")']) 
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dispC ') 

fpriiitf('  Forward  velocity  =  %6.0f  kts\n',Vinf/1.69) 
fprintfC  Temperature  =  %6.0f  degs  F\n',temp) 

fprintfC  Pressure  Altitude  =  %6.0f  ft\n',PA) 
f^rintf('  Gross  weight  =  %6.0f  lbs\n',GW) 

f^rintf('  Number  of  blades  =  %6.0f\n',b) 
fprintfC  Rotor  radius  =  %6.2f  ft\n',R) 

fprintfC  Blade  chord  =  %6.2f  ft\n',cblade) 

fprintfC  Blade  twist  =  %6.2f  degs\n',-l*twist*57.3) 

if  airfoil==l 

dispC  Blade  airfoil  =  FlFl-02') 

elseif  airfoil==2 

dispC  Blade  airfoil  =  VR-12') 

elseif  airfoil==3 

dispC  Blade  airfoil  =  SC1095') 

elseif  airfoil==4 

dispC  Blade  airfoil  =  0012') 

elseif  airfoil==5 

dispC  Blade  airfoil  =  RVR') 

end 

fprintfC  Blade  lift  curve  slope  =  %6.2f  \n',a) 

fprintfC  Blade  wieght  =  %6.2f  lbs\n',wblade) 

fprintfC  Rotational  velocity  =  %6.2frads/sec\n', omega) 

fprintfC  Blade  grip  length  =  %6.2f  ft\n',grip) 

fprintfC  Flinge  offset  =  %6.2f  ft\n',e) 

dispC  ') 

dispC  ') 

dispCpress  any  key  to  continue...') 
pause 

clc 

dispC ') 

dispC  ***  INPUT  DATA  CONTINUED  ***') 
eval(['disp("  ',fdename  1 ,'")']) 

dispC  ') 

fprintfCEquivalent  flat  plate  area  =  %6.2f  fl^2\n',Afh) 
fprintfC  Vertical  Projected  area  =  %6.2f  fU2\n',Afv) 
fprintfC  Wing  Area  =  %6.2f  fU2\n', Swing) 

fprintfC  Wing  Span  =  %6.2f  ft\n',bwing) 

fprintfC  Wing  CL  =  %6.2f  \n',CLwing) 

fprintfC  Wing  CDo  =  %6.4f  \n',CDowing) 

fprintfC  Wing  efficiency  factor  =  %6.2f  \n',ewing) 
fprintfC  Horizontal  tail  area  =  %6.2f  ft^2\n',Shoriz) 

fprintfC  Horizontal  tail  span  =  %6.2f  ft\n',bhoriz) 

fprintfC  Horizontal  tail  CL  =  %6.2f  \n',CLhoriz) 

fprintfC  Horizontal  tail  CDo  =  %6.4f  \n',CDohoriz) 
fprintfC  Vertical  tail  area  =  %6.2f  fU2\n',Svert) 
fprintfC  Vertical  tail  span  =  %6.2f  ft\n',bvert) 
fprintfC  Vertical  tail  CL  =  %6.2f  \n',CLvert) 
fprintfC  Vertical  tail  CDo  =  %6.4f  \n',CDovert) 

fprintfC  Auxiliary  thrust  =  %6.0f  lbs\n',Taux) 
fprintfC  Vertical  Auxiliary  thrust  =  %6.0f  lbs\n',vertaux) 
f^rintfC2/rev  cyclic  input  (deg)  =  %6.0f  \n',theta2c*57.3) 
f^rintfC3/rev  cyclic  input  (deg)  =  %6.2f  \n',theta3s*57.3) 
f^rintf('3/rev  offset  (deg)  =  %6.2f  \n',theta3Poffset*57.3) 
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dispC ') 
disp(' ') 

dispCpress  any  key  to  continue...') 
pause 

clc 

dispC ') 

dispC  ***  CALCULATED  DATA  ***') 
eval(['disp("  ',fdename  1 

dispC  ') 

fprintfC  Fuselage  drag  =  %6.0f  lbs\n',Dfuse) 
fprintfC  Rotor  drag  =  %6.0f  lbs\n',F[rotor) 

l^rintfC  Wing  Lift  =  %6.0f  lbs\n',Lwing) 

l^rintfC  Wing  Drag  =  %6.0f  lbs\n',Dwing) 

fprintfC  Florizontal  tail  lift  =  %6.0f  lbs\n',Lhoriz) 
fprintfC  Florizontal  tail  drag  =  %6.0f  lbs\n',Dhoriz) 
fprintfC  Vertical  tail  side  force  =  %6.0f  lbs\n',Lvert) 
fprintfC  Vertical  tail  drag  =  %6.0f  lbs\n',Dvert) 
l^rintfC  Tip  path  angle  =  %6.2f  degs\n',alphaT*57.3) 
l^rintfC  Rotor  Coning  angle  =  %6.2f  degs\n',betao*57.3) 
fprintfCLocation  of  mean  thmst(r/R)=  %6.2f  \n',rT2) 
l^rintfCCollective  pitch  at  .7  r/R  =  %6.2f  degs\n',thetao*57.3) 
fprintfClst  lat  cyclic  term-Al(deg)=  %6.2f  \n',thetalc*57.3) 
fprintfClst  long  cyclic  term-Bl(deg)=  %6.2f  \n',thetals*57.3) 
dispC  ') 
dispC ') 

dispCpress  any  key  to  continue...') 
pause 

clc 

dispC ') 

dispC  ***  CALCULATED  DATA  CONTINUED  ***') 
eval(['disp("  ',fdename  1 ,'")']) 

dispC  ') 

fprintfC  solidity  (sigma)  =  %6.3f\n', solidity) 

fprintfC  Disk  Loading  =  %6.2f  lbs/fU2\n',DL) 

fprintfC  Figure  of  Merit  =  %6.2f\n',FM) 
l^rintfC  CT/sigma  =  %6.3f  \n',CT_sig) 

l^rintfC  CQ/sigma  =  %6.4f  \n',CQ_sig) 

fprintfC  CH/sigma  =  %6.4f  \n',CFI_sig) 

fprintfCTip  Mach  number  of  adv  blade=  %6.3f  \n',Machtip) 
fprintfC  Advance  ratio  =  %6.3f\n',mu) 

fprintfCRotor  thrust  required  (TPP)=  %6.0f  lbs\n',T) 
fprintfC  Rotor  power  required^  %6.0f  h.p.\n',Protor) 
l^rintfC  Rotor  torque=  %6.0f  ft-lbs\n',Qrotor) 

dispC  ') 
dispC  ') 

dispCpress  any  key  to  continue...') 
pause 

clc 

flag=0; 
else 
dispC ') 

dispC  Enter  a  1  or  2') 
end 
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end 


%  ***  output  to  disk  (text  file)  *** 
dispC ') 

dispC  ***  Saving  data  to  disk  ***') 

pause(2) 

diary  on 

diary  off 

delete  diary 

diary  on 

dispC  ***  RESULTS  ***') 

eval(['disp("  filename  1 

fprintfC  Forward  Velocity  =  %6.0f  kts\n',Vinf/1.69) 

fprintfC  Temperature  =  %6.0f  degs  F\n',temp) 

fprintfC  Pressure  Altitude  =  %6.0f  ft\n',PA) 
l^rintfC  Gross  weight  =  %6.0f  lbs\n',GW) 

l^rintfC  Munber  of  blades  =  %6.0f\n',b) 

fprintfC  Rotor  Radius  =  %6.2f  ft\n',R) 

l^rintfC  Blade  Chord  =  %6.2f  ft\n',cblade) 

l^rintfC  Blade  twist  =  %6.2f  degs\n',-l*twist*57.3) 

l^rintfC  Blade  lift  curve  slope  =  %6.2f  \n',a) 
l^rintfC  Blade  weight  =  %6.2f  lbs\n',wblade) 

fprintfC  Rotational  Velocity  =  %6.2frads/sec\n', omega) 
l^rintfC  Blade  grip  length  =  %6.2f  ft\n',grip) 
fprintfC  Flinge  offset  =  %6.2f  ft\n',e) 
fprintfC  Equivalent  fiat  plate  area=  %6.2f  ft^2\n',Afh) 
fprintfC  Vertical  projected  area  =  %6.2f  fU2\n',AfV) 
fprintfC  Wing  Area  =  %6.2f  fU2\n',Swing) 

fprintfC  Wing  Span  =  %6.2f  ft\n',bwing) 

fprintfC  Wing  CL  =  %6.2f  \n',CLwing) 

fprintfC  Wing  CDo  =  %6.4f  \n',CDowing) 

fprintfC  Wing  Efficicency  Factor  =  %6.2f  \n',ewing) 
fprintfC  Fforizontal  tail  area  =  %6.2f  ft^2\n',Shoriz) 
fprintfC  Fforizontal  tail  span  =  %6.2f  ft\n',bhoriz) 
fprintfC  Fforizontal  tail  CL  =  %6.2f  \n',CLhoriz) 
fprintfC  Fforizontal  tail  CDo  =  %6.4f  \n',CDohoriz) 
fprintfC  Vertical  tail  area  =  %6.2f  fU2\n',Svert) 

fprintfC  Vertical  tail  span  =  %6.2f  ft\n',bvert) 

fprintfC  Vertical  tail  CL  =  %6.2f  \n',CLvert) 

fprintfC  Vertical  tail  CDo  =  %6.4f  \n',CDovert) 
fprintfC  Fuselage  drag  =  %6.0f  lbs\n',Dfuse) 
fprintfC  Rotor  drag  =  %6.0f  lbs\n',Drotor) 

fprintfC  Wing  lift  =  %6.0f  lbs\n',Lwing) 

fprintfC  Wing  drag  =  %6.0f  lbs\n',Dwing) 

fprintfC  Horizontal  tail  lift  =  %6.0f  lbs\n',Lhoriz) 

fprintfC  Horizontal  tail  drag  =  %6.0f  lbs\n',Dhoriz) 
fprintfC  Vertical  tail  side  force=  %6.0f  lbs\n',Lvert) 
l^rintfC  Vertical  tail  drag  =  %6.0f  lbs\n',Dvert) 
fprintfC  Auxiliary  thrust  =  %6.0f  lbs\n',Taux) 
fprintfC  Vertical  Auxiliary  thrust  =  %6.0f  lbs\n',vertaux) 
fprintfC  Tip  path  angle  =  %6.2f  degs\n',alphaT*57.3) 
fj^rintfC  Rotor  coning  angle  =  %6.2f  degs\n',betao*57.3) 
fprintfCLocation  of  mean  thrust(r/R)=  %6.2f  \n',rT2) 
fprintfC  Collective  pitch  at  .7  r/R=  %6.2f  degs\n',thetao*57.3) 
fprintf(Tst  lat  cyclic  term-Al(deg)=  %6.2f  \n',thetalc*57.3) 
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fprintfClst  long  cyclic  term-Bl(deg)=  %6.2f  \n',tlietals*57.3) 
f^rintf('2/rev  cyclic  term  (deg)  =  %6.2f  \n',theta2c*57.3) 
f^rintf('3/rev  cyclic  term  (deg)  =  %6.2f  \n',tlieta3s*57.3) 
f^rintf('3/rev  offset  (deg)  =  %6.2f  \n',tlieta3Poffset*57.3) 
f^rintf('  solidity  =  %6.3f\n', solidity) 

fprintfC  Disk  loadings  %6.2f lbs/ft^2\n',DL) 

fprintfC  Figure  of  Merit  =  %6.2f\n',FM) 

f^rintf('  CT/sigma  =  %6.3f  \n',CT_sig) 

f^rintf('  CQ/sigma  =  %6.4f  \n',CQ_sig) 

f^rintf('  CH/sigma  =  %6.4f  \n',CFt_sig) 

fprintf('Tip  Mach  of  the  adv.  blade  =  %6.3f  \n',Machtip) 
fprintfC  Advance  Ratio  =  %6.3f\n',mu) 

fprintf('Rotor  thrust  required(TPP)  =  %6.0f  lbs\n',T) 
fprintfC  Rotor  power  required  =  %6.0f  h.p.\n',Protor) 
fprintfC  Rotor  Torque  =  %6.0f  ft-lbs\n',Qrotor) 
diary  off 


%  ***  Output  to  disk  (.mat  file  containing  matrix  variables: 

%  r,  psi,  vi,  theta,  thetaXP,  betat,  alpha,  Tpsi,  Mpsi,  DMpsi,  dT,  dM 
%  dD,  CLplot,  CDplot)  *** 

%  ***  Configuring  variables  for  output  *** 

theta=theta*57.3; 

thetaXP=thetaXP*57.3; 

betat=[betat  twist*(0.7-(Reff+(R-Reff)/2)/R)]*57.3; 
alpha=alpham*57.3;,alpha=[alpha  zeros(size(psi))]; 

Mpsi=Mpsi( : ,  length(Mpsi(  1 , : ))- 1 ) ; 

dM=[dM  ddM]; 

psi=psi*57.3; 

r=[r  (R-(R-Reff)/2)]; 

vi=[vi  0]; 

dr=dr*12; 


clc 

flagl=l; 
while  flagl  >  0 
dispC  ') 

wantplots=inputCDo  you  want  to  view  plots,  1 .  yes  or  2.  no?'); 
if  wantplots==2, 
flag  1=0; 

elseif  wantplots==l, 
plotter 
flag  1=0; 
else 
dispC ') 

dispC  Enter  a  1  or  2') 
end 
end 

clc 

dispC ') 

dispC  ***  PERFORMANCE  OUTPUT  DATA  INSTRUCTIONS  ***' 
dispC ') 

dispC  A.  Single  value  data  saved  to  a  file  named:') 
eval(['disp("  '",filenamel,'.prf '")']) 
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dispC ') 

dispC  B.  This  is  a  text  fde,  use  the  tyoe  command  to  view  the') 
dispC  file  or  use  a  text  editor  to  view/print  the  file.') 
eval(['fiag=exist('", filename  l,'.prf');']) 
if  flag  <  1, 

eval(['!rename  diary  ',filenamel,'.prf"]) 
else 

eval([' !  del ', filename  1  ,'.prf "]) 
eval(['!rename  diary  ',filenamel,'.prf"]) 
end 

filename2=[filenamel  '_p']; 
dispC ') 

dispC  C.  Matrix  and  vector  data  saved  to  a  file  named:') 
eval(['disp("  '",filename2,'.mat"")']) 
dispC  ') 

dispC  D.  This  is  a  ".mat"  binary  file,  use  the  load  command') 
dispC  to  retrieve  the  data  for  plotting.') 

eval(['save  ',filename2,'  r  psi  vi  theta  thetaXP  betat  alpha  Tpsi  Mpsi  DMpsi  dT  dM  dD  CLplot 

CDplof]) 

dispC ') 

dispC  ***  END  OF  PERFORMANCE  ***') 
dispC ') 
dispC ') 

dispC  press  any  key  to  continue...') 
pause 
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APPENDIX  D.  PLOTTER.M 


%  ***  Plotting  Routine  added  28  July  2003  some  code  copied/modified  from  lines  28-72  of 

create_plots 

%  JANRAD  v6.0  *** 
clc 

dispC ') 

dispC  ****  Available  Plots  ****') 
dispC ') 

dispC  1 .  Rotor  Azimuth  Angle  v.  Cyclic  Pitch') 
dispC  2.  Blade  Azimuth  Position  v.  Airload') 
dispC  3.  Rotor  Airload  Distribution') 
dispC  4.  Blade  Azimuth  Position  v.  Lift') 
dispC  5.  Rotor  Lift  Distribution') 

dispC  6.  Angle  of  Attack  Distribution') 

dispC  7.  Rotor  Tip  Path  Plane  Stall  pattern  and  Thrust  Distribution') 

dispC  8.  Rotor  Trim  Lissajous  Figure') 

dispC  9.  Lift  Distribution  (Contour  Plot)') 

dispC  10.  CL  Distribution  (Contour  Plot)') 

dispC  1 1 .  Mach  Number  Distribution  (Contour  Plot)') 

dispC ') 

dispC ') 

flag2=l; 

drinches=dr; 

while  flag2  >  0 

plotchoice=input('Please  choose  a  plot  by  entering  1-11.'); 

ifplotchoice==l, 

flag2=0; 

if  twoPinput==0  &  threePinput==0 
figure(l) 
plot(psi,thetaXP) 
axis([0,360,-10,30]); 
legend('lP  Harmonic  Cyclic  Pitch') 

title(['Rotor  Azimuth  Angle  v.  Harmonic  Cyclic  Pitch  for  ',num2str(Vinf/l. 68781),... 

'  Kts,  advance  ratio=',num2str(mu)]) 
xlabel('Rotor  Azimuth  Angle  (degs)'),ylabel('Cychc  Pitch  (degs)') 
elseifrvrCP==l 
figure(l) 

thetao=thetao.*ones(size(theta2Pr)); 

plot(psi,thetao*57.3,'-',psi,(-theta-thetao*57.3),':',psi,theta2Pr*57.3,'-',psi,thetaXP,'*') 

axis([0,360,-25,20]); 

legend('Collective  Pitch', 'IP  Cyclic  Pitch','2P  Cyclic  Pitch','Coll  -i-  IP  -i-  2P  signal') 
title(['Rotor  Azimuth  Angle  v.  RVR  Harmonic  Cyclic  Pitch  for  ',num2str(Vinf/1.68781),... 

'  Kts,  advance  ratio=',num2str(mu),',  2P  input=',num2str(twoPinput),... 

'  degs']) 

xlabel('Rotor  Azimuth  Angle  (degs)'),ylabel('Cychc  Pitch  (degs)') 
else 

figure(l) 

plot(psi,theta,psi,theta2Pr*57.3,'-',psi,theta3Pr*57.3,':',psi,thetaXP,'*') 

axis([0,360,-15,35]); 

legend('lP  Cyclic  Pitch','2P  Cyclic  Pitch',... 

'3P  Cyclic  Pitch','lP  2P  3P  signal') 

title(['Rotor  Azimuth  Angle  v.  Harmonic  Cyclic  Pitch  for  ',num2str(Vinf/1.68781),... 
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'  Kts,  advance  ratio=',num2str(mu),',  2P  input=',num2str(twoPinput),... 

'  degs,  3P  input=',num2str(threePinput),... 

'  degs,  3P  offset=',num2str(threePoffset),'  degs']) 
xlabel('Rotor  Azimuth  Angle  (degs)'),ylabel('Cyclic  Pitch  (degs)') 
end 

elseif  plotchoice==2, 
flag2=0; 
figure(2) 

subplot(3,3,8) 

plot(r./R,dT(l,:)./drinches,'k'),grid 
title('Psi  =  0  deg  ') 

xlabel('Blade  Position  r/R');ylabel('Airload  (Lb/in)'); 
subplot(3,3,2) 

plot(r./R,dT(floor(naz/2),:)./drinches,'k'),grid 
title('Psi  =180  deg ') 

xlabel('Blade  Position  r/R');  ylabel('Airload  (Lb/in)'); 
subplot(3,3,6) 

plot(r./R,dT(floor(naz/4),:)./drinches,'k'),grid 
title('Psi  =  90  deg  ') 

xlabel('Blade  Position  r/R');ylabel('Airload  (Lb/in)'); 
subplot(3,3,4) 

plot(r./R,dT(floor(3*naz/4),:)./drinches,'k'),grid 
title('Psi  =  270  deg  ') 

xlabel('Blade  Position  r/R');  ylabel('Airload  (Lb/in)'); 

gtext(['Blade  Position  v.  Airload  for  ',num2str(Vinf/l. 68781),'  Kts,  advance 

ratio=',num2str(mu),. .. 

',  2P  input=',num2str(twoPinput),'  degs,  3P  input=',... 
num2str(threePinput), '  degs,  3P  offset=',num2str(threePoffset),'  degs']) 

elseif  plotchoice==3, 
flag2=0; 
figure(3) 

[th,rl]  =  meshgrid((-180:360/naz:180)*pi/180,r/R); 

[x,y]  =  pol2cart(th,rl); 
dTl=[dT;  dT(l,:)]; 
for  i=l:naz+l 

dT  1  (i, :  )=dT  1  (i, : )  ./(drinches) ; 
end 

mesh(x',y',dTl) 
view(3 15,60) 

xlabel('Starboard');  ylabel('Aff);  zlabel('Aero  Load,  Lb/in'); 

title(['Airload  Distribution  At  ',num2str(Vinf/l. 68781),'  Kts,  advance  ratio=',num2str(mu),... 

',  2P  input=',num2str(twoPinput),'  degs,  3P  input=',... 
num2str(threePinput), '  degs,  3P  offset=',num2str(threePoffset),'  degs']) 

elseif  plotchoice==4, 
flag2=0; 
figure(4) 
subplot(3,3,8) 
plot(r./R,dN  ( 1 , :  )),grid 
title('Psi  =  0  deg  ') 
xlabel('r/R') ;  ylabel('Lift,Lb') ; 
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subplot(3,3,2) 

plot(r./R,dN(floor(naz/2),:)),grid 
title('Psi  =180  deg ') 
xlabel('r/R') ;  ylabel('Lift,Lb') ; 

subplot(3,3,6) 

plot(r./R,dN(floor(naz/4),:)),grid 
title('Psi  =  90  deg  ') 
xlabel('r/R') ;  ylabel('Lift,Lb') ; 

subplot(3,3,4) 

plot(r./R,dN(floor(3  *naz/4),  :)),grid 
title('Psi  =  270  deg  ') 
xlabel('r/R');  ylabel('Lift,lb'); 

gtext(['Blade  Position  v.  Lift  for',num2str(Vinf/ 1.68781),'  Kts,  advanee  ratio-, num2str(mu),... 
',  2P  input=',num2str(twoPinput),'  degs,  3P  input=',... 
num2str(threePinput), '  degs,  3P  offset=',num2str(threePoffset),'  degs']) 

elseif  plotehoice==5, 
flag2=0; 
figure(5) 

[th,rl]  =  meshgrid((-180:360/naz:180)*pi/180,r/R); 

[x,y]  =  pol2eart(th,rl); 
dNl=[dN;  dN(l,:)]; 
for  i=l:naz+l 
dNl(i,:)=dNl(i,:); 
end 

mesh(x',y',dNl) 
view(3 15,60) 

xlabel('Starboard');  ylabel('Aff);  zlabel('Lift,  Lb'); 

title(['Lifl  Distribution  At  ',num2str(Vinf/L68781),'  Kts,  advanee  ratio- ,num2str(mu),... 

',  2P  input=',num2str(twoPinput),'  degs,  3P  input=',... 
num2str(threePinput), '  degs,  3P  offset=',num2str(threePoffset),'  degs']) 

elseif  plotehoiee==6, 
flag2=0; 
figure(6) 

alpha=alpham*  57.3; 

alpha=alpha'; 

drl=(R-e)/nbe; 

rl=e:drl:R-drl; 

rl=rl+drl/2; 

[th,r2]=meshgrid((0 : 3  60/naz:  3  60)  *pi/ 1 80,r  1  /R) ; 

[x,y]=pol2eart(th,r2); 
alpha=[alpha,alpha(:,  1)]; 

%  The  veetor  v  may  be  altered  to  show  user  defined  eontour  lines 
ifrvrCP==l 

v=[-196,-194,-190,-180,-170,-18,-14,-10,-8,-4,-2,0,2,4,8,10,14,18]; 

else 

v=[18,-16,-14,-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,16]; 

end 

hl=polarl([0,2*pi],  [0,1]); 

delete(hl) 

hold  on 
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hub  area 


%  Either  [C,h2]  equation  may  be  used.  The  second  equation,  shown  in  line  162,  depicts  the 
%  [C,h2]=contour(x,y,alpha,v); 


[C,h2]=contour(x(round(grip*12/dr):nbe,l:nbe+l),y(round(grip*12/dr):nbe,l:nbe+l),alpha(round(grip*12/ 

dr):nbe,l:nbe+l),v); 

clabel(C,h2); 

title(['Angle  of  Attack  Distribution  At  ',num2str(VinEl. 68781),'  Kts,  advance 

ratio=',num2str(mu),. .. 

',  2P  input=',num2str(twoPinput),'  degs,  TPP  angle=',... 

num2str(alphaT*57.3),'  degs,  3P  input=',num2str(threePinput), '  degs,  3P  offset=',... 
num2str(threePoffset),'  degs']) 
hold  on 

thetahub=[0:pi/270:2*pi]; 
hub=(grip/R) ./( 1  +0  *cos(thetahub)) ; 
for  i=l:.  1:20 
polar(thetahub,hub/i); 
hold  on 
end 

hold  off 
view(90,90) 

elseif  plotchoice==7, 
flag2=0; 
figure(7) 

subplot(2,l,l) 

plot(psi,alpha(naz,  l:naz)),set(gca,'XTick', [0:90:360], 'YTick', [-8:2: 14]) 
xlabel('Azimuth  Angle,  degs'),ylabel('Tip  Angle  of  Attack,  degs') 
title(['Rotor  Tip  Path  Plane  Stall  Pattern,  advance  ratio=',num2str(mu)]) 
text(120,10,'Ideal') 
hold  on 

psiplotl=  [0,45,45,135,135,180,210,250,270,290,360]; 
ideal=  [10,7,-.5,-.5,7,10,ll, 12,12.5,13, 10.5]; 
plot(psiplotl,ideal,'r— ') 
hold  off 


dTplot=sum(dT')/(rho*Adisk*(omega*R)^2); 

subplot(2,l,2) 

plot(psi,dTplot),set(gca,'XTick',  [0:90:360]) 

xlabel('Azimuth  Angle,  degs'),ylabel('Thrust  Coefficient  (one  blade)') 
title(['Actual  and  Ideal  Thrust  Distribution  at  Stall,  advance  ratio=',num2str(mu)]) 
text(180,.003,'Ideal') 
hold  on 

psiplot2=  [0  45  45  135  135  180  210  240  270  315  340  360]; 

ideal2=  [.003  .0048  .0001  .0001  .0048  .003  .002  .001  .0009  .0015  .002  .0029]; 
plot(psiplot2,ideal2,'r— ') 
hold  off 

gtext([num2str(Vinf/l. 68781),'  Kts,  advance  ratio=',num2str(mu),... 

',  2P  input=',num2str(twoPinput),'  degs,  TPP  angle=',... 

num2str(alphaT),'  degs,  3P  input=',num2str(threePinput), '  degs,  3P  offset=',... 

num2str(threePoffset),'  degs']) 

elseif  plotchoice==8, 
flag2=0; 
figure(8) 
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Mpsil=[Mpsi;  Mpsi(l,l)]; 

%  if  desired  convert  to  in  *  lbs 
Mpsil=Mpsil*12; 
thetal=[theta;  theta(l,l)]; 
plot(Mpsil,thetal) 
axis([-le6,3e6,0,35]); 
theta  1  =theta  1  *pi/ 1 8  0 ; 
alphaT=alphaT*  5  7 . 3 

%  determine  Lissajous  area 
area=trapz(Mpsi*  12, theta/57. 3); 
area=abs(area); 

xlabel('Thmst  Moment,  in  lb'),ylabel('Blade  Angle,  degs') 

title(['Lissajou  Figure  ',num2str(Vinf/l. 68781),'  Kts,  advance  ratio=',num2str(mu),... 

',  2P  input=',num2str(twoPinput),'  degs,  Lissajou  Area=',num2str(area),'  in  lbs,  TPP 
angle=',num2str(alphaT),'  degs,  3P  input=',... 

num2str(threePinput), '  degs,  3P  offset=',num2str(threePoffset),'  degs']) 

elseif  plotchoice==9, 
flag2=0; 
figure(9) 
drl=(R-e)/nbe; 
rl=e:drl:R-drl; 
rl=rl+drl/2; 

[th,r2]=meshgrid((0 : 3  60/naz:  3  60)  *pi/ 1 80,r  1  /R) ; 

[x,y]=pol2cart(th,r2); 


for  i=l:naz; 
for  j=l:naz; 

dT2(i,j)=dT(i,j)./(drinches); 

end 

end 

dT3=[dT2;  dT2(l,:)]; 
dT3=dT3'; 

v=[-5,-4,-3,-2,-l,0,4,8,12,16,20,24,28,32,36,40,44,48,52,56,60,64,68]; 
hl=polarl([0,2*pi],  [0,1]); 
delete(hl) 
hold  on 

%  Either  [C,h2]  equation  may  be  used.  The  second  depicts  the  hub  area 
%  [C,h2]=contour(x,y,dT3,v); 

[C,h2]=contour(x(round(grip*12/dr):nbe,l:nbe+l),y(round(grip*12/dr):nbe,l:nbe+l),dT3(round(grip*  12/dr 
):nbe,l:nbe+l),v); 

clabel(C,h2); 

title(['Airload  Distribution  At  ',num2str(Vinf/l. 68781),'  Kts,  advance  ratio=',num2str(mu),... 

',  2P  input=',num2str(twoPinput),'  degs,  TPP  angle=',... 

num2str(alphaT*57.3),'  degs,  3P  input=',num2str(threePinput), '  degs,  3P  offset-,... 
num2str(threePoffset),'  degs']) 
hold  on 

thetahub=[0:pi/270:2*pi]; 
hub=(grip/R) ./( 1  +0  *cos(thetahub)) ; 
for  HI:.  1:20 
polar(thetahub,hub/i); 
hold  on 
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end 

hold  off 
view(90,90) 

elseif  plotchoice==10, 
flag2=0; 
figure(lO) 
drl=(R-e)/nbe; 
rl=e:drl:R-drl; 
rl=rl+drl/2; 

[th,r2]=meshgrid((0 : 3  60/naz:  3  60)  *pi/ 1 80,r  1  /R) ; 

[x2,y2]=pol2cart(th,r2); 

CLplotl^CLplof; 

CLplotl=[CLplotl,CLplotl(:,l)]; 

%  The  vector  v  may  be  altered  to  show  user  defined  contour  lines 

%v=[-1,-.8,-.6,-.4,--3,-2,-.1,0,-1,-2,.3,.4,-6,.8,1.0,1.2,1.4]; 

v=[-l, -.8, -.6, -.4, -.2,0, -2, .4, -6, .8, 1.0, 1-2, 1.4]; 

hl=polarl([0,2*pi],  [0,1]); 

delete(hl) 

hold  on 

[C,h2]=contour(x2(round(grip*12/dr):nbe,l:nbe+l),y2(round(grip*12/dr):nbe,l:nbe+l),CLplotl(round(gri 

p*12/dr):nbe,l:nbe+l),v); 

clabel(C,h2); 

title(['CL  Distribution  At  ',num2str(Vinf71. 68781),'  Kts,  advance  ratio=',num2str(mu),... 

',  2P  input=',num2str(twoPinput),'  degs,  TPP  angle=',... 

num2str(alphaT*57.3),'  degs,  3P  input=',num2str(threePinput), '  degs,  3P  offset=',... 
num2str(threePoffset),'  degs']) 
hold  on 

thetahub=[0:pi/270:2*pi]; 
hub=(grip/R) ./( 1  +0  *cos(thetahub)) ; 
for  HI:.  1:20 
polar(thetahub,hub/i); 
hold  on 
end 

hold  off 
view(90,90) 

elseif  plotchoice==l  1, 
flag2=0; 
figure(ll) 
drH(R-e)/nbe; 
rl=e:drl:R-drl; 
rHrl+drl/2; 

[th,r2]=meshgrid((0 : 3  60/naz:  3  60)  *pi/ 1 80,r  1  /R) ; 

[x3  ,y3  ]=pol2cart(th,r2) ; 

machplot  1  =machplot'; 
machplot2=[machplot  1  ,machplot  1  ( : ,  1 )] ; 

v=[-.4,-.3,-.2,-.l,0,.l,.2,.3,.4,.5,.6,.7,.8,.9,l]; 
hl=polarl([0,2*pi],  [0,1]); 
delete(hl) 
hold  on 

[C,h2]=contour(x3,y3,machplot2,v); 
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clabel(C,h2); 

title(['Mach  Number  Distribution  At  ',num2str(Vinf71. 68781),'  Kts, 
ratio=',num2str(mu)]) 
hold  on 

thetahub=[0:pi/270:2*pi]; 
hub=(grip/R) ./( 1  +0  *cos(thetahub)) ; 
for  i=l:.l:20 

polar(thetahub,hub/i); 
hold  on 
end 

hold  off 
view(90,90) 

else 

dispC ') 

dispC  ***  Enter  a  1,  2,  3,  4,  5,  6,  7,  8,  9,  10,  or  11  ***  ') 
end 
end 

clc 

flag3=l; 
while  flag3  >  0 
dispC ') 

wantmoreplots=input('Do  you  want  to  view  additional  plots,  1.  yes  or  2.  no?'); 
if  wantmoreplots==2, 
flag3=0; 

elseif  wantmoreplots==l 
plotter 
flag3=0; 
else 
dispC ') 

dispC  Enter  a  1  or  2') 
end 
end 


%  ***  End  of  Plotting  Routine  *** 


advance 
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APPENDIX  E.  POLARl.M 


function  hpol  =  polar(theta,rho,line_style) 

%POLAR  Polar  coordinate  plot. 

%  POLAR(THETA,  RHO)  makes  a  plot  using  polar  coordinates  of 
%  the  angle  THETA,  in  radians,  versus  the  radius  RHO. 

%  POLAR(THETA,RHO,S)  uses  the  linestyle  specified  in  string  S. 
%  See  PLOT  for  a  description  of  legal  linestyles. 

% 

%  See  also  PLOT,  LOGLOG,  SEMILOGX,  SEMILOGY. 

%  Copyright  1984-2002  The  MathWorks,  Inc. 

%  SRevision:  5.22  $  $Date:  2002/04/08  21:44:28  $ 

if  nargin  <  1 

error('Requires  2  or  3  input  arguments.') 
elseif  nargin  ==  2 
if  isstr(rho) 
linestyle  =  rho; 
rho  =  theta; 

[mr,nr]  =  size(rho); 
ifmr==  1 
theta  =  1  :nr; 
else 

th  =  (l:mr)'; 
theta  =  th(:,ones(l,nr)); 
end 
else 

linestyle  =  'auto'; 
end 

elseif  nargin  ==  1 
linestyle  =  'auto'; 
rho  =  theta; 

[mr,nr]  =  size(rho); 
ifmr==  1 
theta  =  l:nr; 
else 

th  =  (l:mr)'; 
theta  =  th(:,ones(l,nr)); 
end 
end 

if  isstr(theta)  |  isstr(rho) 

error('Input  arguments  must  be  numeric.'); 
end 

if~isequal(size(theta),size(rho)) 

error('THETA  and  RHO  must  be  the  same  size.'); 
end 

%  get  hold  state 
cax  =  newplot; 

next  =  lower(get(cax,'NextPlof)); 
holdstate  =  ishold; 

%  get  x-axis  text  color  so  grid  is  in  same  color 
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tc  =  get(cax,'xcolor'); 

Is  =  get(cax,'gridlinestyle'); 

%  Hold  on  to  current  Text  defaults,  reset  them  to  the 

%  Axes'  font  attributes  so  tick  marks  use  them. 

lAngle  =  get(cax,  'DefaultTextFontAngle'); 

fName  =  get(cax,  'DefaultTextFontName'); 

fSize  =  get(cax,  'DefaultTextFontSize'); 

fWeight  =  get(cax,  'DefaultTextFontWeight'); 

fUnits  =  get(cax,  'DefaultTextUnits'); 

set(cax,  'DefaultTextFontAngle',  get(cax,  'FontAngle'),  ... 

'DefaultTextFontName',  get(cax,  'FontName'),  ... 

'DefaultTextFontSize',  get(cax,  'FontSize'), ... 

'DefaultTextFontWeight',  get(cax,  'FontWeight'),  ... 

'DefaultT  extUnits',  'data') 

%  only  do  grids  if  hold  is  off 
if  -holdstate 

%  make  a  radial  grid 
hold  on; 

maxrho  =  max(abs(rho(:))); 

hhh=plot([-maxrho  -maxrho  maxrho  maxrho], [-maxrho  maxrho  maxrho  -maxrho]); 
set(gca,'dataaspectratio',[l  1  l],'plotboxaspectratiomode','auto') 

V  =  [get(cax,'xhm')  get(cax,'ylim')]; 
ticks  =  sum(get(cax,'ytick')>=0); 
delete(hhh); 

%  check  radial  limits  and  ticks 

rmin  =  0;  rmax  =  v(4);  rticks  =  max(ticks-l,2); 
if  rticks  >  5  %  see  if  we  can  reduce  the  number 
if  rem(rticks,2)  ==  0 
rticks  =  rticks/2; 
elseifrem(rticks,3)  ==  0 
rticks  =  rticks/3; 
end 
end 

%  define  a  circle 
th  =  0:pi/50:2*pi; 
xunit  =  cos(th); 
yunit  =  sin(th); 

%  now  really  force  points  on  x/y  axes  to  he  on  them  exactly 
inds  =  l:(length(th)-l)/4:length(th); 
xunit(inds(2:2:4))  =  zeros(2,l); 
yunit(inds(l:2:5))  =  zeros(3,l); 

%  plot  background  if  necessary 
if  ~isstr(get(cax,'color')), 
patch('xdata',xunit*rmax,'ydata',yunit*rmax,  ... 
'edgecolor',tc,'facecolor',get(gca,'color'),... 

'handlevisibihty','off); 

end 

%  draw  radial  circles 

c82  =  cos(82*pi/180); 
s82  =  sin(82*pi/180); 
rticks=l; 
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rinc  =  (rmax-rmiii)/rticks; 

for  i=(rmin+riiic):rinc:rmax 
hhh  =  plot(xuiiit*i,yuiiit*i, Is, 'color', tc,'linewidth',l,... 

'handlevisibility',  'off) ; 
text((i+rinc/20) *c82, (i+rinc/20) * s82,  ... 

['  '  num2str(i)],'verticalalignment','bottom',... 
'bandlevisibility','off) 
end 

set(bbb,'bnestyle','-')  %  Make  outer  circle  solid 

%  plot  spokes 

tb  =  (l:6)*2*pi/12; 

cst  =  cos(tb);  snt  =  sin(tb); 

cs  =  [-cst;  cst]; 

sn  =  [-snt;  snt]; 

plot(rmax*cs,rmax*sn,ls,'color',tc,'bnewidtb',l,... 

'bandlevisibibty','off) 

%  annotate  spokes  in  degrees 
rt  =  1 . 1  *rmax; 
for  i  =  1  :lengtb(tb) 

text(rt*cst(i),rt*snt(i),int2str(i*30),... 

'borizontalalignmenf ,  'center', . . . 

'bandlevisibibty','off); 
if  i  ==  lengtb(tb) 
loc  =  int2str(0); 
else 

loc  =  int2str(180-l-i*30); 
end 

text(-rt*cst(i),-rt*snt(i),loc,'borizontalalignment','center',... 

'bandlevisibibty','off) 

end 

%  set  view  to  2-D 
view(2); 

%  set  axis  limits 

axis(rmax*[-l  1  -1.15  1.15]); 
end 

%  Reset  defaults. 

set(cax,  'DefaultTextFontAngle',  fAngle  ,  ... 
'DefaultTextFontName',  fName  ,  ... 

'DefaultTextFontSize',  fSize, ... 

'DefaultTextFontWeigbf,  fWeigbt,  ... 
'DefaultTextUnits',fUnits ); 

%  bansform  data  to  Cartesian  coordinates. 

XX  =  rbo.*cos(tbeta); 
yy  =  rbo.*sin(tbeta); 

%  plot  data  on  top  of  grid 
if  strcmp(line_style,'auto') 
q  =  plot(xx,yy); 
else 

q  =  plot(xx,yy,line_style); 
end 
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if  nargout  >  0 
hpol  =  q; 
end 

if  -holdstate 

set(gca,'dataaspectratio',[l  1  1]),  axis  off;  set(cax,'NextPlof,next); 
end 

set(get(gca,'xlabel'),'visibleVon') 
set(get(gca,'ylabel'), 'visible', 'on') 
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APPENDIX  F.  THRCALC.M 


%  thrcalc  calculates  the  total  thrust  along  a  blade  at 
%  each  azimuth  (psi)  location 
%  dispCbegin  thrcalc') 

Up=zeros(size(psi  *r)) ; 

Ut=zeros(size(Up)); 

dT=zeros(size(Up)); 

ddT=zeros(size(psi)); 

dN=zeros(size(Up)); 

ddN=zeros(size(psi)); 

for  i=l  :length(psi), 

Up(i,:)=vi.*cos(betao)+Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 
Ut(i,:)=r.*omega+Vinf*cos(alphaT)*sin(psi(i)); 
phi=atan2(Up(i, :  ),Ut(i, : )) ; 

mach  =  ( Vtip .  *cos(alphaT) .  *r./R+ V inf  *  sin(psi(i)))/(49 . 05  *  sqrt(temp+460)) ; 

%  alpha=theta(i)+betat-phi; 
alpha=thetaXP(i)+betat-phi; 


ifairfoil==l, 

[CL,CD]=hh02clcd(alpha); 

%  additional  airfoils  added  for  experimentation 
elseif  airfoil==2, 

[CL,CD]=vrl2clcd(alpha); 
elseif  airfoil==3, 

[CL,CD]=sc  1 095r8clcd(alpha,mach); 
elseif  airfoil==4, 

[CL,CD]=oo  1 2clcd(alpha,mach); 
elseif  airfoil==5, 

[CL,CD]=rvrclcd(alpha,mach); 

end 

dT(i,:)=0.5*rho.*cblade.*dr.*(Up(i,:).^2+Ut(i,:).^2).*(CL.*cos(phi)-CD.*sin(phi)); 
Tpsi(i)=sum(dT  (i, :)); 

dN(i,:)=0.5*rho.*cblade.*dr.*(Up(i,:).^2+Ut(i,:).^2).*(CL.*cos(alpha)+CD.*sin(alpha)); 
Npsi(i)=sum(dT  (i, : )) ; 

%  Populate  CLplot  and  CDplot  matrix 
CLplot(i,:)=CL; 

CDplot(i,:)=CD; 

%  Populate  Mach  matrix 
machplot(i, :  )=mach; 

%  Populate  alpha  matrix 
alpham(i,  :)=alpha; 

%  ***  calculations  for  tip  loss  area  *** 

Uptip=Vinf*sin(alphaT)*cos(betao)+VinPcos(alphaT)*sin(betao)*cos(psi(i)); 

Uttip=(R-(R-Reff)/2)*omega+Vinf*cos(alphaT)*sin(psi(i)); 

phitip=atan2(Uptip,Uttip); 

alphatip=theta(i)+betat-phitip; 

%ddT(i)=0.5*rho*cblade*(R-Reff)*(Uptip''2)*(-.009*sin(phitip)); 
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ddT(i)=0.5*rho*cblade*(0.5+0.5*cos(2*psi(i)))*(R-Reff)*(Uptip^2+Uttip^2)*(- 

.009*siii(phitip)); 

T  psi(i)=Tpsi(i)+ddT  (i) ; 

ddN(i)=0.5*rho*cblade*(0.5+0.5*cos(2*psi(i)))*(R-Reff)*(Uptip^2+Uttip''2)*(- 

.009*siii(alphatip(i))); 

Npsi(i)=Npsi(i)+ddN  (i) ; 
end 

macbplot2=macbplot'; 

CLplot2=CLplot'; 
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APPENDIX  G.  TMCALC.M 


%  tmcalc  calculates  the  total  thrust  moment  along  a  blade 
%  at  each  azimuth  (psi)  location 
%  dispCbegin  tmcalc') 

Up=zeros(size(psi  *r)) ; 

Ut=zeros(size(Up)); 

dM=zeros(size(Up)); 

ddM=zeros(size(psi)); 

for  i=l  :length(psi), 

Up(i,:)=vi.*cos(betao)+Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Ut(i,:)=r.*omega+Vinf*cos(alphaT)*sin(psi(i)); 

phi=atan2(Up(i,:),Ut(i,:)); 

%  alpha=theta(i,k)+betat-phi; 
alpha=thetaXP(i,k)+betat-phi; 
ifairfoil==l, 

[CL,CD]=hh02clcd(alpha); 

%  additional  airfoils  added  for  experimentation 
elseif  airfoil==2, 

[CL,CD]=vrl2clcd(alpha); 
elseif  airfoil==3, 

[CL,CD]=sc  1 095r8clcd(alpha,mach); 
elseif  airfoil==4, 

[CL,CD]=oo  1 2clcd(alpha,mach); 
elseif  airfoil==5, 

[CL,CD]=rvrclcd(alpha,mach); 

end 

%  Populate  CLplot  and  CDplot  matrix 
CLplot(i,:)=CL; 

CDplot(i,:)=CD; 

%  Populate  alpham  matrix 
alpham(i,  :)=alpha; 


dM(i,:)=0.5*rho*cblade.*r*dr.*(Up(i,:).^2+Ut(i,:).^2).*(CL.*cos(phi)-CD.*sin(phi)); 
Mpsi(i,k)=sum(dM(i, : )) ; 

%  ***  calculations  for  tip  loss  areas  *** 

Uptip=Vinf*sin(alphaT)*cos(betao)+VinPcos(alphaT)*sin(betao)*cos(psi(i)); 

Uttip=(R-(R-Reff)/2)*omega+Vinf*cos(alphaT)*sin(psi(i)); 

phitip=atan2(Uptip,Uttip); 

ddM(i)=0.5*rho*cblade*(R-(R-Reff)/2)*(R-Reff)*(Uptip^2+Uttip^2)*(-.009*sin(phitip)); 

Mpsi(i,k)=Mpsi(i,k)+ddM(i); 

end 

CLplot2=CLplof; 
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APPENDIX  H.  TRIM.M 


%  Rotor  trim  subroutine 
dispC  ') 

dispC  ***  EXECUTING  ROTOR  TRIM  ROUTINE  ***') 


%  ***  calculation  of  required  parameters  *** 

rlio=.002377*(-.00003  l*PA+(-.002*temp+l .  118)); 

%  Hardwire  rho  for  high  and  hot  conditions 
if  PA==4000  &  temp==95 
rho=.00192 
else 
end 


%  ***  first  guess  at  rotor  profile  drag  (  H  force  )  *** 

ifVinf<  16.9, 

Drotor=0; 

else 

Drotor=VinP(rho/.002377); 

end 

q=0 . 5  *  rho  *  VinU2 ; 

Adisk=pi*R*2; 

Vtip=omega*R; 

Dfuse=q*Afh; 

CDwing=CDowing+(CLwing^2/(ewing*pi*(bwing^2/Swing))); 
CDhoriz=CDohoriz+(CLhoriz^2/(.  8  *pi  *  (bhoriz^2/Shoriz))) ; 

CDvert=CDo  vert+(CLvert''2/(.  8  *pi  *  (bvert^2/S  vert))) ; 

Dwing=q*CDwing*Swing; 

Dhoriz=q*CDhoriz*  Shoriz; 

Dvert=q*CDvert*  Svert; 

Lwing=q*CLwing*  Swing; 

Lhoriz=q*CLhoriz*  Shoriz; 

Lvert=q*CLvert*  Svert; 

LftotaULwing+Lhoriz; 

DflotaUDfuse+Dwing+Dhoriz+Dvert; 

%  ***  thrust  calculation  *** 

%  Check  percentage  of  wing  lift  versus  rotor  thrust 

if  helochoice==2  &  Vinf  >  16.9 
clc 

LoverT=LftotaEG  W ; 

disp('The  wing  is  carrying  '),disp(LoverT*100),disp('percent  of  the  aircraft  weight.') 
flag=l; 

while  flag  >  0 
disp(") 

disp('Do  you  want  to  change  the  percentage  of  aircraft  weight  carried  by  the  wing?  ') 
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disp('l.  continue  with  current  wing  settings  or  2.  change  percentage  of  aircraft  weight 
carried  by  the  wing  by  varying  wing  area') 

wingchoice=input('Enter  a  1  or  2...'); 
while  isempty(wingchoice), 
dispC  ') 

disp('You  must  enter  a  1  or  2') 

disp('Do  you  want  to  change  the  percentage  of  aircraft  weight  carried  by  the  wing?  ') 
disp('l.  continue  with  current  wing  settings  or  2.  change  percentage  of  aircraft  weight 
carried  by  the  wing  by  varying  wing  area') 

wingchoice=input('Enter  a  1  or  2...'); 
end 

ifwingchoice~=l  &  wingchoice~=2 
dispC  ') 

disp('  ***  Enter  a  1  or  2  ***') 
elseif  wingchoice==2 

wpercent=input('Enter  desired  percentage  of  aircraft  weight  carried  by  wing:  '); 

lifterror=wpercent/100*GW-Lflotal; 

if  lifterror  >  0 

while  lifterror  >  .05*wpercent/100*GW 
Swing=Swing  +  5; 

Lwing=q*CLwing*  Swing; 

Lftotal=Lwing+Lhoriz; 

lifterror=wpercent/100*GW-Lftotal; 

end 

lifterror=0; 

disp('New  wing  area  is  '),disp(Swing),disp('  ft^2') 

disp('Press  any  key  to  continue  ...') 

pause 

elseif  lifterror  <  0 

while  lifterror  <  .05*wpercent/100*GW 
Swing=Swing  -  5; 

Lwing=q*CLwing*  Swing; 

Lftotal=Lwing+Lhoriz; 

lifterror=wpercent/100*GW-Lftotal; 

end 

lifterror=0; 

disp('New  wing  area  is  '),disp(Swing),disp('  ft^2') 
disp('Press  any  key  to  continue  ...') 
pause 
end 
flag=0; 

elseif  wingchoice==l 

dispC  Wing  settings  unchanged') 
disp('Press  any  key  to  continue  ...') 
pause 
flag=0; 
end 
end 


end 


%  if  required,  define  alphaT  fixed  quantity 
if  alphaTfix==l 

alphaT=alphaTinput; 

else 

alphaT =atan2((Dftotal+Drotor),  (GW-Lftotal)); 
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end 


%  Determine  thrust  parameters  for  low  speed  flight 
ifVinf<  16.9, 

%  Determine  augmented  vertical  thrust  capability 
clc 

disp('JANRAD  is  about  to  perform  low  speed  performance  analysis.') 
disp(") 
flag=l; 
while  flag  >  0 
disp(") 

disp('Are  any  auxiliary  vertical  thmst  devices  installed  on  the  aircraft,') 

disp('i.e..  Lift  Fans  or  Direct  Jet  thrusters?') 

disp(") 

vertauxchoice=input('l.  yes  or  2.  no'); 
while  isempty(vertauxchoice), 
dispC  ') 

disp('You  must  enter  a  1  or  2') 

disp('Are  there  any  auxiliary  vertical  thmst  devices  installed  on  the  aircraft,') 

disp('i.e..  Lift  Fans  or  Direct  Jet  thmsters?') 

dispC) 

vertauxchoice=input('l.  yes  or  2.  no'); 
end 

ifvertauxchoice~=l  &  vertauxchoice~=2 
dispC  ') 

dispC  ***  Enter  a  1  or  2  ***') 
elseif  vertauxchoice==l 

vertaux=inputCEnter  amount  of  thmst  in  pounds  provided  by  the  vertical  auxiliary  thmst 

devices: '); 

GW=GW  -vertaux; 

dispCInput  parameters  adjusted  to  reflect  auxiliary  vertical  thmst.') 
flag=0; 

elseif  vertauxchoice==2 
vertaux=0; 

dispCNo  vertical  thmst  devices  installed') 
dispCPress  any  key  to  continue  ...') 
pause 
flag=0; 
end 
end 

T=(l+(0.3*Afv/Adisk))*GW; 

CT=T/(Adisk*rho*Vtip^2); 

%  Determine  thmst  parmeters  for  nominal  flight  speeds 
else 

T=(GW-Lftotal)/cos(alphaT); 

CT=T/(Adisk*rho*Vtip^2); 

vertaux=0; 


%  Added  for  Compound/RVR  auxiliary  thmst  analysis 
if  helochoice==2  &  Vinf  >  16.9 
clc 

if  T  auxchoice==  1 

T  aux=(Dfuse+Dwing+Dhoriz+Dvert) ; 

%  if  required,  define  alphaT  fixed  quantity 

113 


if  alphaTfix==l 

alphaT=alphaT  input ; 
else 

alphaT=atan2((Dftotal+Drotor),  (GW-Lftotal)); 
end 
else 

if  Taux<((Dfuse+Dwing+Dhoriz+Dvert)-sin(alpliaT)*T) 

dispC***  Auxiliary  thrust  will  not  overeome  drag  at  this  airspeed,  ***') 
dispC  Total  drag  =  '),disp(Dfuse+Dwing+Dhoriz+Dvert) 

dispC  Thrust  defieit  ='),disp(Dfuse+Dwing+Dhoriz+Dvert-sin(alphaT)*T-Taux) 
flag=l; 
while  flag  >  0 
disp(") 

disp('Do  you  want  to  ehange  the  auxiliary  thrust  value?  ') 

disp(T.  eontinue  with  eurrent  auxiliary  thrust  value  or  2.  set  auxiliary  thrust  to  equal 

drag') 

Tauxehoiee=input('Enter  a  1  or  2...'); 
while  isempty(Tauxehoiee), 
dispC  ') 

dispC  You  must  enter  a  1  or  2') 
disp(") 

dispCDo  you  want  to  ehange  the  auxiliary  thrust  value?  ') 
disp(") 

dispC  1.  eontinue  with  eurrent  auxiliary  thrust  value  or  2.  set  auxiliary  thrust  to  equal 

drag') 

Tauxehoiee=inputCEnter  a  1  or  2...'); 
end 

if  Tauxehoiee~=l  &  Tauxehoiee~=2 
dispC ') 

dispC  ***  Enter  a  1  or  2  ***') 
elseif  Tauxehoiee==2 

T  aux=(Dfuse+Dwing+Dhoriz+Dvert) ; 
flag=0; 

dispCAuxillary  thrust  ehange  to  equal  total  drag') 
disp(") 

dispCPress  any  key  to  eontinue  ...') 
pause 

elseif  T  auxehoiee==  1 

dispCAuxillary  thrust  defieit  is  '),disp(Dfuse+Dwing+Dhoriz+Dvert-sin(alphaT)*T- 

Taux), dispC  lbs.') 

dispC***  Results  may  be  erroneous.  ***') 
disp(") 

dispCPress  any  key  to  eontinue  ...') 
pause 
flag=0; 
end 
end 
end 

%  define  alphaT  fixed  quantity 
if  alphaT fix==l 

alphaT=alphaT  input ; 
else 

alphaT =atan2((Dftotal+Drotor),  (GW-Lftotal)); 
end 
end 
end 
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Dftotal=(Dfuse+Dwing+Dhoriz+Dvert)-Taux; 

mu=ViiiPcos(alphaT)A^tip; 

optmu-0107*(Vinf/1.69)-1.2; 


%  if  desired,  iterate  omega  to  obtain  desired  mu  value 
if  heloehoiee==2 
ifrvrCP— 1 
ifmu<  1.0 
ele 

dispCThe  advanee  ratio,  mu,  is'),disp(mu) 
disp(") 
flag=l; 
while  flag  >  0 

disp('This  value  is  too  low  for  RVR  operation,  do  you  want  JANRAD  to  optimize  mu  so 

the') 

disp('rotor  system  will  operate  in  RVR  mode?') 

disp('l.  eontinue  with  eurrent  mu  value  or  2.  ehange  mu  by  automatieally  varying 
rotational  speed'); 

muehoiee=input('Enter  a  1  or  2...'); 
while  isempty(muehoiee), 
dispC ') 

disp('You  must  enter  a  1  or  2') 

disp('Do  you  want  JANRAD  to  optimize  mu  so  the  rotor  system  will  operate  in  RVR 

mode? ') 

disp('l.  eontinue  with  eurrent  mu  value  or  2.  ehange  mu  by  automatieally  varying 

rotational  speed'); 

muehoiee=input('Enter  a  1  or  2...'); 
disp('Optimum  mu  value  eurrently  set  to'),disp(optmu) 
end 

if  muehoiee~=l  &  muehoiee~=2 
dispC ') 

dispC  ***  Enter  a  1  or  2  ***') 
elseif  muehoiee==2 
muerroi=mu-optmu; 
if  muerror  <  0 

while  abs(muerror)  >  .01 
omega=omega  -  .001; 

Vtip=omega*R; 
mu=V  inPeos(alphaT)  tip ; 
muerroi^mu-optmu; 
end 

muerror=0; 

disp('New  rotational  velocity  is  '),disp(omega) 
disp(") 

disp('New  mu  value  is  '),disp(mu) 
disp(") 

disp('Press  any  key  to  continue  ...') 
pause 
else 
end 
flag=0; 

elseif  muchoice==l 

disp('Rotational  velocity  and  mu  value  unchanged.') 
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disp(") 

disp('Press  any  key  to  continue  ...') 
pause 
flag^O; 
end 
end 
else 
end 
else 
end 
else 
end 

ifmu<=  1.0&rvrCP==l 
clc 

dispCThe  advance  ratio  is  less  than  1.0  and  you  have  chosen  RVR  based') 
disp('cychc  pitch.  JANRAD  cannot  complete  computations  using  this  input') 
disp('The  program  will  revert  to  conventional  Higher  Harmonic  Feathering.') 
disp('Press  any  key  to  continue  ...') 
pause 
rvrCP=2; 
else 
end 

ifrvrCP==l  &  twoPinput==0 
clc 

disp('You  have  selected  RVR  cyclic  pitch,  2/rev  cyclic  input  cannot  be  zero.') 
twoPinput=input('Please  enter  a  non-zero  2/rev  cyclic  input  in  degrees  or  press  <CONTROL  -l- 
C>  to  quit  JANRAD'); 

while  isempty(twoPinput), 
dispC ') 

disp('You  must  enter  a  value') 
twoPinput=input('2/rev  cyclic  input  (degs): '); 
while  twoPinput==0 

disp('The  value  cannot  be  zero  either!') 

twoPinput=input('Please  enter  a  non-zero  2/rev  cyclic  input  or  press  <CONTROL  -I-  C>  to 

quit  JANRAD'); 

end 

end 

else 

end 

clc 

dispC  ***  ANALYSIS  IN  PROGRESS  ***') 

%  ***  setup  blade  radius  elements,  azimuth  elements, 

%  induced  velocity  distributions,  and  determination 
%  of  coning  angle  and  tip  loss  parameters  *** 

B=l-(sqrt(2*CT)/b); 

Reff=B*R; 

Rbar=Reff-e; 
dr=(Reff-grip)/nbe; 
r=grip:dr:Reff-dr;,r=r-l-dr/2; 
rTl=.7;,%  ***  first  guess  at  rT  *** 

RbarT=rT  1  *Rbar; 
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mblade=wblade/32. 1 7 ; 

betao=asin((T/b*RbarT-(.5*(R-e)+e)*wblade)/((.5*(R-e)+e)^2*omega''2*mblade)); 

betat=twist*(0.7-(r/R)); 

psi=0:360/naz:360-360/ naz ;  ,p  si=psi  '/5  7 . 3 ; 

%  ***  induced  velocity  distribution  *** 

ifVinf<  16.9, 

A=4*pi; 

Bv=(b/2)*omega*a*cblade; 

Tv=0; 

delT=T-Tv; 

while  abs(delT)  >  .01*T 

thetav=twist*  (0 . 7-(r/R))+thetao; 

C=(-b/2)*cblade*omega^2*r*a.*thetav; 

vi=(-Bv+sqrt(Bv^2-(4*A*C)))/(2*A); 

dTv=(b/2)*rho*((omega*r).^2)*a.*(thetav-(vi./(omega*r)))*cblade*dr; 

Tv=sum(dTv); 

delT=T-Tv; 

ifdelT<0, 

thetao=thetao-0.5*thetao*abs(delT/T); 

else 

thetao=thetao+0.5*thetao*abs(delT/T); 

end 

end 

else 

lamdaT=[  1  -2*mu*sin(alphaT)  mu^2*(sin(alphaT)^2+l)-2*mu^3  *sin(alphaT) 

mu^4*sin(alphaT)^2-CT^2/4]; 

lamdaT  =max(real(roots(lamdaT))) ; 
vi=lamdaT  *Vtip-Vinf*sin(alphaT); 
vi=vi*ones(size(r)); 
end 

%  ***  first  guess  at  theta  *** 

thetalc=0.035*((0.0006e-3*Vinf^2+0.244e-3*Vinf)/0.105); 
thetals=-0.087*((0.0006e-3*Vinf^2+0.244e-3*Vinf)/0.105); 
theta2c=twoPinput/57.3 ; 
theta3  s=threePinput/5  7.3; 
theta3Poffset=threePoffset/57.3 ; 

theta=thetao+theta  1  c .  *cos(psi)+theta  Is.*  sin(psi) ; 

theta2Pr=theta2c.  *cos(2*psi); 

theta3  Pr=theta3  s.  *  sin(3  *psi+theta3  Poffset) ; 

thetaXP=thetao+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta 

3Poffset); 

ifrvrCP==l 

thetaXP=thetao+thetao*sin(psi)+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)+theta2c.*cos(2*psi); 

else 

end 


%  ***  rotor  trimming  routine  *** 

dispC  ') 
disp(' ') 
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dispC  ***  TRIMMING  COLLECTIVE  ***') 
k=l; 


ifrvrCP— 1 
devTC=.20; 
devTC- 1; 
else 

devTC=.02; 

end 

errorO=(T*devTC)+ 1 ; 
while  abs(errorO)  >  T*devTC; 

Tpsi=zeros(size(psi)); 

threale 

en'orO=T-(mean(Tpsi)  *b) ; 
if  errorO  <  -T*devTC 

thetao=thetao-0.35*thetao*abs(L5*error0/T)*abs(l-mu);%  added  abs  for  RVR  operation 
elseif  errorO  >  T*devTC 

thetao=thetao+0.35*thetao*abs(L5*error0/T)*abs(l-mu);%  added  abs  for  RVR  operation 
end 

theta=thetao+theta  1  e .  *eos(psi)+theta  Is.*  sin(psi) ; 

theta2Pr=theta2e.  *eos(2*psi); 

theta3  Pr=theta3  s.  *  sin(3  *psi+theta3  Poffset) ; 

thetaXP=thetao+thetale.*eos(psi)+thetals.*sin(psi)+theta2e.*eos(2*psi)+theta3s.*sin(3*psi+theta3Poffset) 

ifrvrCP==l 

thetaXP=thetao+thetao*sin(psi)+thetale.*eos(psi)+thetals.*sin(psi)+theta2e.*eos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetale.*eos(psi)-thetals.*sin(psi)+theta2e.*eos(2*psi); 

else 

end 

ifk>  1, 

if  abs(errorO)  >=  abs(errorl), 
ele 

dispC ') 
dispC ') 

disp('This  eonfiguration  will  not  trim') 
disp('Try  a  lower  airspeed  or  a  new  design') 
dispC  ') 

disp('Program  exeeution  terminated') 
dispC  ') 

errorC***  END  OF  PROGRAM  ***') 
end 
end 

error  l=errorO; 
ifrvrCP— 1 

errorO=(T*devTC)-l ; 
end 
end 

dispC ') 

dispC  ***  TRIMMING  CYCLIC  ***') 
tO=eloek; 
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k=l; 


if  rvrCP==2 

e^•orO=(((T^)  *rT  1  *  (R-grip))  * .  04)+ 1 ; 
while  errorO  >  ((T/b)*rTl*(R-grip))*.04 
time=etime(clock,tO) ; 

if  time  >15, 
dispC ') 

disp('still  trimming') 
t0=cloek; 
end 

Mpsi(:,k)=zeros(size(psi)); 

tmeale 

theta=[theta  theta(:,k)]; 
theta2Pr=[theta2Pr  theta2Pr( :  ,k)] ; 
theta3  Pr=[theta3  Pr  theta3  Pr( :  ,k)] ; 
thetaXP=[thetaXP  thetaXP(:,k)]; 

Mpsi=[Mpsi  Mpsi(:,k)]; 

%  ***  ealeulation  of  initial  dthetadM  *** 

ifk  <  2, 

theta( :  ,k+ 1  )=theta( :  ,k)+0 . 25/5 7. 3 ; 
theta2Pr( :  ,k+ 1  )=theta2Pr( :  ,k)+0 .25/57.3; 
thetaXP(:,k+l)=thetaXP(:,k)+0.25/57.3; 

Mpsi( :  ,k+ 1  )=zeros(size(psi)) ; 

k=k+l; 

tmeale 

k=k-l; 

dthetadM=(thetaXP( :  ,k+ 1  )-thetaXP( :  ,k))  ./(Mpsi( :  ,k+ 1  )-Mpsi( :  ,k)) ; 
end 

%  ***  ealeulation  of  M  first  harmonie  parameters  *** 

M 1  e=2  *  sum(Mpsi( :  ,k) .  *eos(psi))/naz; 

M 1  s=2  *  sum(Mpsi  (:,k).*sin(psi))/ naz ; 

%  ***  removal  of  first  harmonie  terms  from  Mpsi  *** 

Mpsi(:  ,k+ 1  )=Mpsi( :  ,k)-M  1  e .  *eos(psi)-M  1  s .  *  sin(psi) ; 
delM=Mpsi( :  ,k+ 1  )-Mpsi( :  ,k) ; 
errorO=max(delM)-min(delM); 
ifk>  1, 

if  errorO  >  error  1, 
ele 

dispC ') 
dispC ') 

disp('This  eonfiguration  will  not  trim') 
disp('Try  a  lower  airspeed  or  a  new  design') 
disp(' ') 

disp('Program  exeeution  terminated') 
disp(' ') 

errorC***  END  OF  PROGRAM  ***') 
end 
end 
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error  l=errorO; 


%  ***  calculation  of  new  theta  *** 

delM=0 . 5  *  ( 1  -mu)  *delM; 
theta( :  ,k-i- 1  )=theta( :  ,k)-i-(dthetadM .  *  delM) ; 
theta2Pr( :  ,k-l- 1  )=theta2Pr( :  ,k)-l-(dthetadM.  *delM) ; 
theta3  Pr( :  ,k-i- 1  )=theta3  Pr( :  ,k)-i-(dthetadM .  *  delM) ; 
thetaXP( :  ,k-i- 1  )=thetaXP( :  ,k)-i-(dthetadM.  *delM) ; 

if  errorO  <=  ((T/b)*rTl*(R-grip))*.04, 
theta  1  c=2  *  sum(theta( :  ,k) .  *cos(psi))/naz; 
theta  1  s=2  *  sum(theta( :  ,k) .  *  sin(psi))/naz; 

%  ***  2P  coefficient  *** 

theta2c=sum(ones(size(theta( :  ,k))) .  *twoPinput/ 57.3  )/naz; 

%  ***  3P  coeeficient  *** 

theta3  s=sum(ones(size(theta( :  ,k))) .  *threePinput/ 57.3  )/naz; 
else 

theta  1  c=2*  sum(theta( :  ,k-l- 1 ) .  *cos(psi))/naz; 
theta  1  s=2  *  sum(theta( :  ,k-l- 1 ) .  *  sin(psi))/naz; 

%  ***  2P  coefficient  *** 

theta2c=sum(ones(size(theta( :  ,k-l- 1 ))) .  *twoPinput/ 57.3  )/naz; 

%  ***  3P  coefficient  *** 

theta3  s=sum(ones(size(theta( :  ,k-l- 1 ))) .  *threePinput/5  7. 3  )/naz; 
end 

%  ***  calculate  theta  with  IP  terms  only  *** 

theta( :  ,k-i- 1  )=thetao-i-theta  1  c .  *cos(psi)-i-theta  Is.*  sin(psi) ; 

%  ***  add  2P  and  3P  terms,  result  thetaXP  *** 
theta2Pr( :  ,k-l- 1  )=theta2c  .*cos(2*psi); 
theta3  Pr( :  ,k-l- 1  )=theta3  s .  *  sin(3  *psi-l-theta3  Poffset) ; 

thetaXP(:,k-l-l)=thetao-l-thetalc.*cos(psi)-l-thetals.*sin(psi)-l-theta2c.*cos(2*psi)-t-theta3s.*sin(3*psi-l-theta3 

Poffset); 

ifrvrCP— 1 


thetaXP(:,k-l-l)=thetao-t-thetao*sin(psi)-l-thetalc.*cos(psi)-l-thetals.*sin(psi)-l-theta2c.*cos(2*psi); 

thetaXP=thetao-l-thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)-l-theta2c.*cos(2*psi); 

else 

end 

%  ***  calculation  of  new  dthetadM  *** 
theta=[theta  theta(:,k-l-l)]; 
theta2Pr=[theta2Pr  theta2Pr( :  ,k-l- 1 )] ; 
theta3  Pr=[theta3  Pr  theta3  Pr( :  ,k-l- 1 )] ; 
thetaXP=[thetaXPthetaXP(:,k-H)]; 

Mpsi=[Mpsi  Mpsi(:,k-l-l)]; 

theta( :  ,k-i-2)=theta( :  ,k)-i-0 . 25/5 7. 3 ; 
theta2Pr(:,k-K2)=theta2Pr(:,k)-K0.25/57.3; 
theta3Pr(:,k-K2)=theta3Pr(:,k)-K0.25/57.3; 
thetaXP(:,k-K2)=thetaXP(:,k)-K0.25/57.3; 

Mpsi(:  ,k-i-2)=zeros(size(Mpsi(:  ,k-i- 1 ))); 
k=k-l-2; 
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tmcalc 

k=k-2; 

dthetadM=(thetaXP( :  ,k+2)-thetaXP( :  ,k)) ./ (Mpsi( :  ,k+2)-Mpsi( :  ,k)) ; 
k=k+l; 
end 
else 

Mpsi(:,k)=zeros(size(psi)); 

tmcalc 

theta=[theta  theta(:,k)]; 
tlieta2Pr=[tlieta2Pr  tlieta2Pr(:  ,k)] ; 
thetaS  Pr=  [thetaS  Pr  thetaS  Pr( :  ,k)] ; 
thetaXP=[thetaXP  thetaXP(:,k)]; 

Mpsi=[Mpsi  Mpsi(:,k)]; 

%  *  *  *  calculation  of  initial  dthetadM  *  *  * 

ifk  <  2, 

tlieta( :  ,k+ 1  )=theta( :  ,k)+0 .25/5 7. 3 ; 

theta2Pr( :  ,k+ 1  )=theta2Pr( :  ,k)+0 .25/57 . 3 ; 

tlietaXP(:,k+l)=tlietaXP(:,k)+0.25/57.3; 

Mpsi(:,k+l)=zeros(size(psi)); 

k=k+l; 

tmcalc 

k=k-l; 

dtlietadM=(tlietaXP( :  ,k+ 1  )-tlietaXP( :  ,k))  ./(Mpsi( :  ,k+ 1  )-Mpsi( :  ,k)) ; 
end 
end 

dispC ') 

dispC  ***  ADJUSTING  COLLECTIVE  ***') 
disp(' ') 

tlieta=tlieta(:,k); 
theta2Pr=theta2Pr( :  ,k) ; 
theta3  Pr=theta3  Pr( :  ,k) ; 
thetaXP=tlietaXP  ( :  ,k) ; 

k=l; 


ifrvrCP— 1 
devAC=.l; 
devAC=.03; 
else 

devAC=.01; 

end 

errorO=(T*devAC)+l ; 
while  abs(errorO)  >  T*devAC 
Tpsi=zeros(size(psi)); 
thrcalc 

errorO=T-(mean(Tpsi)  *b) ; 
if  errorO  <  -T*devAC, 

thetao=thetao-0.25*thetao*abs(1.25*error0/T)*abs(l-mu); 
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elseif  errorO  >  T*devAC, 

thetao=thetao+.025*thetao*abs(1.25*error0/T)*abs(l-mu); 

end 

tbeta=tbetao+tbeta  1  c.  *cos(psi)+tbeta  Is.*  sin(psi) ; 

tbeta2Pr=tbeta2c .  *cos(2  *psi) ; 

tbeta3  Pr=tbeta3  s .  *  sin(3  *psi+tbreePoffset) ; 

tbetaXP=tbetao+tbetalc.*cos(psi)+tbetals.*sin(psi)+tbeta2c.*cos(2*psi)+tbeta3s.*sin(3*psi+tbeta3Poffset) 


ifrvrCP==l 

tbetaXP=tbetao+tbetao*sin(psi)+tbetalc.*cos(psi)+tbetals.*sin(psi)+tbeta2c.*cos(2*psi); 

tbetaXP=tbetao+tbetao*sin(psi)-tbetalc.*cos(psi)-tbetals.*sin(psi)+tbeta2c.*cos(2*psi); 

else 

end 

ifk>  1, 

if  abs(errorO)  >  abs(errorl), 
clc 

dispC  ') 
dispC  ') 

disp('Tbis  configuration  will  not  trim') 
disp('Try  a  lower  airspeed  or  a  new  design') 
dispC  ') 

disp('Program  execution  terminated') 
dispC  ') 

errorCEND  OF  PROGRAM') 
end 
end 

error  l=errorO; 
k=k+l; 
ifrvrCP— 1 
errorO=(T*devAC)- 1 ; 
end 
end 


%  ***  calculating  drag  moments  *** 

DMpsi=zeros(size(psi)); 

dmcalc 

%  ***  calculating  rotor  FI  force  *** 

ifVinf<  16.9, 

Flrotor=0; 
dT=[dT  ddT]; 
dN=[dN  ddN]; 
dD=[dD  ddD]; 
else 

dT=[dT  ddT]; 
dN=[dN  ddN]; 
dD=[dD  ddD]; 
for  i=l  :lengtb(r)+l, 

Fllc(i)=2*sum(dT(:,i).*cos(psi))/naz; 

FI  1  s(i)=2  *  sum(dT  ( :  ,i) .  *  sin(psi))/naz; 
end 

if  rvrCP==2 

Flrotor=(((b*cos(alpbaT)/2)*(sum(Flls)-sin(betao)*sum(Fllc)))+Drotor)/2; 
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else 

Hrotor=0; 

end 

end 


%  *  *  *  calculating  new  rT  *  *  * 

rT2=(((mean(Mpsi(:  ,lengtli(Mpsi(  1 , 1  ))/mean(Tpsi))/R)+rT  1  )/2; 
%  ***  check  rotor  drag  and  rT,  retrim  rotor  if  needed  *** 


while  abs(Drotor-Hrotor)  >  0.2*Hrotor  |  abs(rTl-rT2)  >  0.5*rTl  %  increase  tolerance  to  50% 
if  abs(Drotor-Hrotor)  >  0.2*Hrotor, 
dispC ') 

dispC  ***  ADJUSTING  ROTOR  DRAG  ***') 
end 


ifrvrCP— 1 
Drotor=0; 
else 

Drotor=Hrotor; 

end 


ifrvrCP— 1 
devMT=.25; 
devMT=.045; 
else 

devMT-015; 

end 

if  abs(rTl-rT2)  >  devMT*rTl, 
dispC ') 

dispC  ***  ADJUSTING  MEAN  THRUST  LOCATION  ***') 
end 

dispC ') 

dispC  ***  RETRIMMING  ROTOR  ***') 
dispC ') 

dT=dT(:,I:nbe); 

dN=dN(:,I:nbe); 

dD=dD(:,l:nbe); 

%  ***  recealculating  parameters  *** 

%  if  required,  define  alphaT  fixed  quantity 
if  alphaTfix==I 

alphaT=aIphaTinput; 

else 

alphaT=atan((Dftotal+Drotor)/(GW-Lftotal)); 

end 

mu=VinPcos(alphaT)A^tip; 

if  Vinf  >=  16.9, 

T=(GW-Lftotal)/cos(alphaT); 


123 


CT=T/(Adisk*rho  *  Vtip^2) ; 

lamdaT=[  1  -2*mu*siii(alphaT)  mu^2*(sin(alphaT)*2+l)-2*mu^3  *siii(alphaT) 

mu^4*siii(alphaT)^2-CT^2/4]; 

lamdaT=max(real(roots(lamdaT))); 

vi=lamdaT*Vtip-VinPsiii(alphaT); 

vi=vi*ones(size(r)); 

end 

B=l-(sqrt(2*CT)/b); 

Reff^B*R; 

Rbar=Reff-e; 
dr=(Reff-grip)/nbe; 
r=grip  :dr :  Reff-dr ;  ,r=r+dr/2 ; 

RbarT=rT2*Rbar; 

betao=asin((T/b*RbarT-(.5*(R-e)+e)*wblade)/((.5*(R-e)+e)''2*omega''2*mblade)); 

%  ***  trimming  collective  *** 

tO=clock; 

k=l; 

errorO=(T*devTC)+l ; 

while  abs(erTorO)  >  T*devTC 
Tpsi=zeros(size(psi)); 
thrcalc 

errorO=T-(mean(Tpsi)*b); 
if  errorO  <  -T*devTC, 

thetao=thetao-0 . 3  5  *thetao  *abs(  1 . 5  *  errorO/T)  *abs(  1  -mu) ; 
elseif  errorO  >  T*devTC, 

thetao=thetao-l-0.35*thetao*abs(1.5*error0/T)*abs(l-mu); 

end 

theta=thetao-l-theta  1  c.  *cos(psi)-l-theta  Is.*  sin(psi) ; 

theta2Pr=theta2c .  *cos(2  *psi) ; 

theta3  Pr=theta3  s .  *  sin(3  *psi-l-threePoffset) ; 

thetaXP=thetao-l-thetalc.*cos(psi)-l-thetals.*sin(psi)-l-theta2c.*cos(2*psi)-f-theta3s.*sin(3*psi-l-theta3Poffset) 


ifrvrCP==l 

thetaXP=thetao-l-thetao*sin(psi)-l-thetalc.*cos(psi)-l-thetals.*sin(psi)-f-theta2c.*cos(2*psi); 

thetaXP=thetao-l-thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)-l-theta2c.*cos(2*psi); 

else 

end 

ifk>  1, 

if  abs(errorO)  >  abs(errorl), 
clc 

dispC  ') 
dispC  ') 

disp('This  configuration  will  not  trim') 
disp('Try  a  lower  airspeed  or  a  new  design') 
dispC  ') 

disp('Program  execution  terminated  at  line  746') 
dispC  ') 

errorC***  END  OF  PROGRAM  ***') 
end 
end 

error  l=errorO; 
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k=k+l; 
ifrvrCP— 1 

errorO=(T*devTC)- 1 ; 
end 
end 

%  ***  trimming  eye  lie  *** 


ifrvrCP— 1 
devCY=.2; 
devCY=.12; 
else 

devCY=.04; 

end 

k=l; 

errorO=(((T  ^)  *rT2  *  (R-grip))  *devC  Y)+ 1 ; 


if  rvrCP==2 

while  errorO  >  ((T/b)*rT2*(R-grip))*devCY 
time=etime(elock,tO); 
if  time  >15, 

dispCstill  trimming...') 
dispC  ') 
tO=eloek; 
end 

Mpsi(:,k)=zeros(size(psi)); 

tmeale 

theta=[theta  theta(:,k)]; 
theta2Pr=[theta2Pr  theta2Pr( :  ,k)] ; 
theta3  Pr=[theta3  Pr  theta3  Pr( :  ,k)] ; 
thetaXP=[thetaXP  thetaXP(:,k)]; 

Mpsi=[Mpsi  Mpsi(:,k)]; 

%  ***  ealeulation  of  initial  dthetadM  *** 

if  k  <  2, 

theta(:  ,k+ 1  )=theta( :  ,k)+0 .25/5 7 .3 ; 

theta2Pr( :  ,k+ 1  )=theta2Pr( :  ,k)+0 .25/5 7 . 3 ; 

theta3Pr( :  ,k+ 1  )=theta3  Pr( :  ,k)+0 .25/5 7 . 3 ; 

thetaXP(:,k+l)=thetaXP(:,k)+0.25/57.3; 

Mpsi(:,k+l)=zeros(size(psi)); 

k=k+l; 

tmeale 

k=k-l; 

dthetadM=(thetaXP( :  ,k+ 1  )-thetaXP( :  ,k))  ./(Mpsi( :  ,k+ 1  )-Mpsi( :  ,k)) ; 
end 

%  ***  ealeulation  of  M  first  harmonie  parameters  *** 

M 1  c=2  *sum(Mpsi( :  ,k) .  *eos(psi))/naz; 

M 1  s=2  *  sum(Mpsi( :  ,k) .  *  sin(psi))/ naz; 

%  ***  removal  of  first  harmonie  terms  from  Mpsi  *** 

Mpsi(:,k+l)=Mpsi(:,k)-Mle.*eos(psi)-Mls.*sin(psi); 
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delM=Mpsi( :  ,k+ 1  )-Mpsi( :  ,k) ; 
errorO=max(delM)-min(delM) ; 
ifk>  1, 

if  errorO  >  error  1, 
clc 

dispC  ') 

disp('This  configuration  will  not  trim') 
disp('Try  a  lower  airspeed  or  a  new  design') 
dispC  ') 

disp('Program  Execution  Terminated  at  line  814') 
disp(' ') 

errorC***  END  OF  PROGRAM  ***') 
end 
end 

error  l=errorO; 

%  ***  calculation  of  new  theta  *** 

delM=0.5*(l-mu)*delM; 
theta( :  ,k+ 1  )=theta( :  ,k)+(dthetadM.  *delM) ; 
theta2Pr( :  ,k+ 1  )=theta2Pr( :  ,k)+(dthetadM.  *delM) ; 
theta3  Pr( :  ,k+ 1  )=theta3  Pr( :  ,k)+(dthetadM.  *delM) ; 
thetaXP( :  ,k+ 1  )=thetaXP( :  ,k)+(dthetadM.  *delM) ; 

if  errorO  <=  ((T/b)*rT2*(R-grip))*.04, 
theta  1  c=2  *  sum(theta( :  ,k) .  *cos(psi))/naz; 
theta  1  s=2  *  sum(theta( :  ,k) .  *  sin(psi))/naz; 

%  2P/3P  coefficients 

theta2c=sum(ones(size(theta( :  ,k))) .  *twoPinput/ 57.3  )/naz; 
theta3  s=sum(ones(size(theta( :  ,k))) .  *  threePinput/ 5  7 .3  )/naz; 
else 

theta  1  c=2*  sum(theta( :  ,k+ 1 ).  *cos(psi))/naz; 
theta  1  s=2  *  sum(theta( :  ,k+ 1 ) .  *  sin(psi))/naz; 

%  2P/3P  coefficients 

theta2c=sum(ones(size(theta( :  ,k+ 1 ))) .  *twoPinput/ 57.3  )/naz; 
theta3  s=sum(ones(size(theta( :  ,k+ 1 ))) .  *  threePinput/57 . 3  )/naz; 
end 

%  calculate  theta  with  IP  terms  only 

theta( :  ,k+ 1  )=thetao+theta  1  c .  *cos(psi)+theta  Is.*  sin(psi) ; 

%  add  2P  and  3P  terms,  unit  value  only  for  theta2c,  result  thetaXP 
theta2Pr( :  ,k+ 1  )=theta2c.  *cos(2  *psi) ; 
theta3Pr(:,k+ 1  )=theta3  s.  *sin(3  *psi+threePoffset); 

thetaXP(:,k+l)=thetao+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta3 

Poffset); 


if  rvrCP==l 


thetaXP(:,k+l)=thetao+thetao*sin(psi)+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP(:,k+l)=thetao+thetao*sin(psi)-thetalc.*cos(psi)- 
theta  Is.*  sin(psi)+theta2c .  *cos(2  *psi) ; 
else 
end 

%  *  *  *  calculation  of  new  dthetadM  *  *  * 
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theta=[theta  theta(:,k+l)]; 
theta2Pr=[theta2Pr  theta2Pr( :  ,k+ 1 )] ; 
theta3  Pr=[theta3  Pr  theta3Pr( :  ,k+ 1 )] ; 
thetaXP=[thetaXP  thetaXP(:,k+l)]; 

Mpsi=[Mpsi  Mpsi(:,k+1)]; 
theta( :  ,k+2)=zeros(size(Mpsi( :  ,k+ 1 ))) ; 
theta2Pr( :  ,k+2)=zeros(size(Mpsi( :  ,k+ 1 ))) ; 
theta3  Pr( :  ,k+2)=zeros(size(Mpsi( :  ,k+ 1 ))) ; 
thetaXP( :  ,k+2)=zeros(size(Mpsi( :  ,k+ 1 ))) ; 

k=k+2; 

tmcalc 

k=k-2; 

%  replace  IP  dthetadM  with  2P  dthetadM 
%  dthetadM=(theta( :  ,k+2)-theta( :  ,k)) ./ (Mpsi( :  ,k+2)-Mpsi( :  ,k)) ; 

dthetadM=(thetaXP( :  ,k+2)-thetaXP( :  ,k))  ./(Mpsi( :  ,k+2)-Mpsi( :  ,k)) ; 
k=k+l; 
end 

else 

Mpsi(:,k)=zeros(size(psi)); 

tmcalc 

theta=[theta  theta(:,k)]; 
theta2Pr=[theta2Pr  theta2Pr( :  ,k)] ; 
theta3  Pr=[theta3  Pr  theta3  Pr( :  ,k)] ; 
thetaXP=[thetaXP  thetaXP(:,k)]; 

Mpsi=[Mpsi  Mpsi(:,k)]; 

%  ***  calculation  of  initial  dthetadM  *** 

ifk  <  2, 

theta( :  ,k+ 1  )=theta( :  ,k)+0 . 25/5 7. 3 ; 
theta2Pr(:  ,k+ 1  )=theta2Pr( :  ,k)+0 .25/57.3; 
theta3  Pr(:  ,k+ 1  )=theta3  Pr( :  ,k)+0 .25/57.3; 
thetaXP(:,k+l)=thetaXP(:,k)+0.25/57.3; 

Mpsi( :  ,k+ 1  )=zeros(size(psi)) ; 

k=k+l; 

tmcalc 

k=k-l; 

dthetadM=(thetaXP(:,k+l)-thetaXP(:,k))./(Mpsi(:,k+l)-Mpsi(:,k)); 

end 

end 

%  ***  retrimming  collective  *** 

theta=theta(:,k); 
theta2Pr=theta2Pr( :  ,k) ; 
theta3  Pr=theta3  Pr( :  ,k) ; 
thetaXP=thetaXP( :  ,k) ; 
k=l; 

errorO=(T*devAC)+l ; 
while  abs(errorO)  >  T*devAC 
Tpsi=zeros(size(psi)); 
thrcalc 
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errorO=T-(mean(Tpsi)*b); 
if  errorO  <  -T*devAC, 

thetao=thetao-0.25*thetao*abs(1.25*error0/T)*abs(l-mu); 
elseif  errorO  >  T*devAC, 

thetao=thetao+0.25*thetao*abs(1.25*error0/T)*abs(l-mu); 

end 

theta=thetao+theta  1  c.  *cos(psi)+theta  Is.*  sin(psi) ; 

theta2Pr=theta2c .  *cos(2  *psi) ; 

theta3  Pr=theta3  s .  *  sin(3  *psi+threePoffset) ; 

thetaXP=thetao+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta3Poffset) 


ifrvrCP==l 

thetaXP=thetao+thetao*sin(psi)+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)+theta2c.*cos(2*psi); 

else 

end 

ifk>  1, 

if  abs(errorO)  >  abs(errorl), 
clc 

dispC  ') 

disp('This  configuration  will  not  trim') 
disp('Try  a  lower  airspeed  or  a  nex  design') 
dispC  ') 

error('END  OF  PROGRAM  at  line  933') 
end 
end 

error  l=errorO; 
k=k+l; 
ifrvrCP— 1 
errorO=(T*devAC)- 1 ; 
end 
end 

%  ***  recalculating  rotor  FI  force  *** 

ifVinf<  16.9, 

Hrotor=0; 
dT=[dT  ddT]; 
dN=[dN  ddN]; 
dD=[dD  ddD]; 
else 

dT=[dT  ddT]; 
dN=[dN  ddN]; 
dD=[dD  ddD]; 
for  i=l  :length(r)+l, 

FI  1  c(i)=2  *  sum(dT  ( :  ,i) .  *cos(psi))/naz; 

FI  1  s(i)=2  *  sum(dD( :  ,i) .  *  sin(psi))/naz; 
end 


ifrvrCP— 1 
Flrotor=0; 
else 

Flrotor=(((b*cos(alphaT)/2)*(sum(Flls)-sin(betao)*sum(Fllc)))+Drotor)/2; 

end 
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end 


%  *  *  *  recalculating  rT  *  *  * 
rTl=rT2; 

rT2=(((mean(Mpsi( :  ,lengtli(Mpsi(  1 , : 1  ))/mean(Tpsi))/R)+rT  1  )/2 ; 
end 


%  ***  recalculating  drag  moments  *** 

dT=dT(:,l:nbe); 

dN=dN(:,l:nbe); 

dD=dD5,l:nbe); 

DMpsi=zeros(size(psi)); 

dmcalc 

dT=[dT  ddT]; 

dN=[dN  ddN]; 

dD=[dD  ddD]; 

clc 

dispC ') 

dispC  ***  ROTOR  TRIMMED  ***') 
dispC  press  any  key  to  continue...') 

pause 
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APPENDIX  I.  0012CLCD.M 


function  [CL,CD]=ool2clcd(alpha,  mach) 

% 

%  ool2clcd  calculates  CL  and  CD  for  the  NACA  0012 
%  airfoil  given  angle  of  attack  in  radians  and  the 
%  local  Mach  number: 

% 

%  [CL,CD]=ool2clcd(alpha,  mach) 

% 

%  Both  'alpha'  and  'mach'  are  intended  to  be  vectors 
%  the  elements  of  which  correspond  to  the  rotor  blade 
%  radial  stations  of  interest  in  a  blade  element  analysis. 
%  All  equations  are  based  on  Ray  Prouty's  treatment  of 
%  the  0012  in  his  text. 


CL=zeros(size(alpha)); 
CD=zeros(size(alpha)); 
a=alpha*  180/pi; 
mach=abs(mach) ; 
aL  =  15  -  16.*mach; 
aD  =  17  -  23.4.*mach; 

K1  =  0.0233  +  0.342.*(mach.^7.15); 
K2  =  2.05  -  0.95.*mach; 


%  CL  for  Mach  numbers  <  0.725  and  AOA  inside  +/-  20  deg: 

chk=(mach<0.725  &  a>=0  &  a<=aL); 

CL=CL+chk.*((0.1./sqrt(l-mach.^2)  -  0.01.*mach).*a); 

chk=(mach<0.725  &  a>aL  &  a<=20); 

CL=CL+chk.*((0.1./sqrt(l-mach.^2)  -  0.01.*mach).*a  -  Kl.*(a-aL).^K2); 
chk=(mach<0.725  &  a>=-20  &  a<-aL); 

CL=CL-chk.*((0.1./sqrt(l-mach.^2)  -  0.01.*mach).*abs(a)  -  Kl.*(abs(a)-aL).^K2); 

chk=(mach<0.725  &  a>=-aL  &  a<0); 

CL=CL-chk.*((0.1./sqrt(l-mach.^2)  -  0.01.*mach).*abs(a)); 


%  CL  for  Mach  numbers  >  0.725  and  AOA  inside  +/-  20  deg: 


chk=(mach>=0.725  &  a>=0  &  a<=aL); 

CL=CL+chk.*((0.677  -  0.744.*mach).*a); 

chk=(mach>=0.725  &  a>aL  &  a<=20); 

CL=CL+chk.*((0.677  -  0.744.*mach).*a  -  (0.0575-0. 144.*(mach-0.725).^0.44).*(a-aL).^(K2)); 

chk=(mach>=0.725  &  a<0  &  a>=-aL); 

CL=CL-chk.*((0.677  -  0.744. *mach).*abs(a)); 
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chk=(mach>=0.725  &  a<-aL  &  a>=-20); 

CL=CL-chk.*((0.677  -  0.744.*mach).*abs(a)  -  (0.0575-0. 144.*(mach-0.725).^0.44).*(abs(a)- 
aL).^(K2)); 


%  CL  for  all  Macb  numbers  and  AOA  outside  +!-  20deg: 

ebk=(a>20  &  a<=161); 

CL=CL-l-ebk.  *(1.15.*  sin(2 .  *alpba)) ; 

ebk=(a>161  &a<=173); 

CL=CL-Kebk.*(-0.7); 

ebk=(a>173  &a<=180); 

CL=CL-Kebk.*(0. 1  .*(a- 1 80)); 


ebk=(a>=-180&a<=-173); 

CL=CL-Kebk.*(0.1.*(a-H80)); 


ebk=(a>-173  &a<=-161); 
CL=CL-Kebk.*(0.7); 

ebk=(a>-161  &  a<-20); 

CL=CL-l-ebk.  *(1.15.*sin(2 .  *alpba)) ; 


%  CD  for  Maeb  numbers  <  0.725  and  AOA  inside  +!-  20  deg: 
ebk=(maeb<0.725  &  a>=0  &  a<=aD); 

CD=CD-Kebk.*(0.0081  (-350.*a  396.*a.^2  -  63.3.*a.^3  3.66.*a.M).*10.''(-6)); 

ebk=(maeb<0.725  &  a>aD  &  a<=20); 

CD=CD-Kebk.*((0.0081  (-350.*a  396.*a.^2  -  63.3.*a.^3  3.66.*a.M).*10.''(-6)) 

0.00066.*(a-aD).^2.54); 

ebk=(maeb<0.725  &  a<0  &  a>=-aD); 

CD=CD-tebk.*(0.0081  (-350.*abs(a)  396.*a.^2  -  63.3.*abs(a).^3  3.66.*a.M).*10.''(-6)); 

ebk=(maeb<0.725  &  a<-aD  &  a>=-20); 

CD=CD-Kebk.*((0.0081  (-350.*abs(a)  396.*a.^2  -  63.3.*abs(a).^3  3.66.*a.M).*10.^(-6)) 

0.00066.*(abs(a)-aD).^2.54); 


%  CD  for  Maeb  numbers  >  0.725  and  AOA  inside  +!-  20  deg: 


ebk=(maeb>=0.725  &  a>=0  &  a<=20); 

CD=CD-Kebk.*((0.0081  (-350.*a  396.*a.^2  -  63.3.*a.^3  3.66.*a.M).*10.''(-6)) 

0.00035. *a.^2.54  21.*(maeb-0.725).^3.2); 

ebk=(maeb>=0.725  &  a<0  &  a>=-20); 

CD=CD-Kebk.*((0.0081  (-350.*abs(a)  396.*a.^2  -  63.3.*abs(a).^3  3.66.*a.M).*10.^(-6)) 

0.00035. *abs(a).^2.54  21.*(maeb-0.725).^3.2); 
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%  CD  for  all  Mach  numbers  and  AOA  outside  +/-  20deg: 


chk=(a>20  &  a<=180); 

CD=CD+chk.*(1.03  -  1.02.*cos(2.*alpha)); 


chk=(a>=-180  &  a<-20); 
CD=CD+chk.*(1.03  -  1.02.*cos(2.*alpha)); 
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APPENDIX!.  HH02CLCD.M 


function  [CL,CD]=hh02clcd(alpha) 

%  hh02clcd  calculates  CL  and  CD  for  an  HH-02  airfoil 
%  given  angle  of  attack  (alpha)  in  radians  and  Mach  number  (Mach) 
%  [CL,CD]=hh02clcd(alpha,mach) 

CL=zeros(size(alpha)); 

CD=zeros(size(alpha)); 
a=alpha*  180/pi; 


chkl=(a>=20  &a<=180); 

CL=CL+chkl.*(0.42541+0.026863*a+5.5988e-4*a.^2-2.1493e-5*a.^3+1.5932e-7*a.M-3.4659e 

10*a.^5); 

CD=CD+chkl.*(-0.7179+0.061213*a-5.9861e-4*a.^2+7.3708e-6*a.^3-6.6605e-8*a.M+1.913e- 

10*a.^5); 


chkl=(a>=-180  &  a<=-50); 

CL=CL+chkl. *(-4.6183-0. 1923*a-3.5554e-3*a.^2-3.3273e-5*a.^3-1.4528e-7*a.M-2.3003e- 
10*a.^5); 

CD=CD-Kchkl.*(2.7093e-2-2.1309e-2*a-K2.0335e-4*a.^2-K3.47e-7*a.^3-3.0586e-8*a.M-1.2584e 

10*a.^5); 

chkl=(a>-50  &  a<-20); 

CL=CL-Kchkl.*(-2.5519-0.22847*a-9.5667e-3*a.^2-1.7051e-4*a.^3-1.0909e-6*a.M); 

CD=CD-Kchkl.*(2.7093e-2-2.1309e-2*a-K2.0335e-4*a.^2-K3.47e-7*a.^3-3.0586e-8*a.M-1.2584e 

10*a.^5); 


chkl=(a>=-20  &  a<=-10); 

CL=CL-Kchkl.*(-0.2-K0.089*a-K0.0034*a.^2); 

CD=CD-Kchkl.*(2.7093e-2-2.1309e-2*a-K2.0335e-4*a.^2-K3.47e-7*a.^3-3.0586e-8*a.M-1.2584e 

10*a.^5); 

chkl=(a<20  &  a>-10); 

CL=CL-Kchkl.*(5.8766e-2-H.3131e-l*a-K2.4742e-3*a.^2-5.303e-4*a.^3-1.5818e-5*a.M-M.28e- 
6*a.^5); 
chk2=a<-4; 
chk2=chk2.*chkl ; 

CD=CD-Kchk2.*(1.3786-^0.916*a-K0.21396*a.^2-K2.0371e-2*a.^3-^7.0076e-4*a.M); 
chk2=(a>=-4  &  a<=7); 
chk2=chk2.*chkl ; 

CD=CD-Kchk2.*(9.732e-3-^3.2326e-4*a-M.4392e-4*a.^2-8.5073e-5*a.^3-H.1826e- 
6*a.M-H.5271e-6*a.^5); 
chk2=a>7; 
chk2=chk2.*chkl ; 

CD=CD-Kchk2.*(1.842e-l-5.7532e-2*a-K5.8043e-3*a.^2-1.2803e-4*a.^3); 
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APPENDIX  K.  RVRCLCD.M 


function  [CL,CD]=rvrclcd(alpha,mach) 

%RVRCLCD  calculates  CL  and  CD  for  a  12%  thickness  Fairchild 
%  RVR  airfoil  given  angle  of  attack  (alpha)  in  radians  and 
%  mach  number  (mach).  Reference  Fairchild  Report  FtC144R1070 
%  figs  D-28  and  D-12 

%  JANRAD  ver  lP2P3Pbeta 

CL=zeros(size(alpha)); 

CD=zeros(size(alpha)); 

Mach_CL=[0  .26  .4  .6  .8  .9]; 

aoa_CL=[-180  -150  -120  -90  -30  -24  -20  -16  -12  -8  -4  0  4  8  12  16  20  24  30  90  120  150  160  170 
180]; 


0  0 

0 

0 

0 

0; 

0 

.7 

.7 

.7  .7 

.7 

•> 

0 

.4 

.4 

.4  .4 

.4 

•> 

0 

0 

0 

0  0 

0; 

0 

-.75 

-.75 

-.75 

-.75 

-.75; 

0 

-.89 

-.89 

-.89 

-.89 

-.89; 

0 

-.9 

-.9 

-.9  -. 

9  -, 

.9; 

0 

-1.0 

-1.0 

-1.0 

-1.0 

-1.0; 

0 

-1.05 

-1.05  -1.05  -1.05  -1.05; 

0 

-.85 

-.85 

-.85 

-.85 

-.85; 

0 

-.5 

-.5 

-.5  -. 

5  -, 

.5; 

0 

.25 

.25 

.25 

.25 

.25; 

0 

.25 

.25 

.25 

.25 

.25; 

0 

.7 

.7 

.7  .5 

.55; 

0 

.9 

.9 

.9  .7 

.65; 

0 

1.52 

1.52 

:  1.52 

.78 

1.0; 

0 

1.65 

1.65 

:  1.65 

.7 

.95; 

0 

1.5 

1.5 

1.5 

1.2 

1.2; 

0 

1.7 

1.7 

1.7 

1.7 

1.7; 

0 

0 

0 

0  0 

0; 

0 

-.55 

-.55 

-.55 

-.55 

-.55; 

0 

-1.09 

-1.09  -1.09  -1.09  -1.09; 

0 

-.8 

-.8 

-.8  -. 

65  - 

-.6; 

0 

-1.1 

-1.1 

-1.1 

-.7 

-.5; 

0 

0 

0 

0  0 

0] 

•> 

Mach_CD=[0  .26  .4  .6  .8  .9]; 

aoa_CD=[-180  -150  -120  -90  -30  -24  -20  -16  -12  -8  -4  0  4  8  12  16  20  24  30  90  120  150  160  170 


180]; 
cd  tab=[ 

0  0 

0 
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.2 

.22; 

0 

.4 

.4 

.4 

.4 

.4; 

0 

.42 

.42 

.42 

.42 

.42; 

0 

.5 

.5 

.5 

.5 

.5; 

0 

.6 

.6 

.6 

.6 

.6; 

0 

1.65 

1.65  1.65  1. 

,65  1. 

0 

1.2 

1.2 

1.2 

1.2 

1.2; 

0 

.65 

.65 

.65 

.65 

.65; 

0 

.45 

.45 

.45 

.45 

.45; 

0 

.15 

.15 

.17 

.2 

.21; 

0 

0 

0 

0 

.05 

.10]; 

rvra=alpha.  *  1 80/pi; 
forj  =  l:length(rvra) 
ifrvra(j)  <  -180 
rvra(j)  =  rvra(j)  -1-360; 
elseif  rvra(j)  >180 
rvra(j)  =  rvra(j)  -  360; 
end 
end 

maeh  =  abs(maeh); 

CL  =  interp2(Maeh_CL,aoa_CL,el_tab,maeh,rvra); 
CD  =  interp2(Maeh_CD,aoa_CD,ed_tab,maeh,rvra); 


138 


APPENDIX  L.  SC1095R8CLCD.M 


function  [CL,CD]=sc  1 095r8clcd(alpha,mach) 

%scl095r8  calculates  CL  and  CD  for  a  Sikorsky  SC1095R8  airfoil 
%given  angle  of  attack  (alpha)  in  radians  and  mach  number  (mach) 
%[CL,CD]=scl095r8clcd(alptia,mach) 


CL=zeros(size(alpha)); 

CD=zeros(size(alpha)); 


Mach_CL=[0  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0]; 

aoa_CL=[-180  -178  -176  -174  -172  -170  -168  -166  -164  -162  -160  -158  -90  -22  -20  -18  -16  -14  - 
12  -10  -8  ... 

-6  -4  -2  0  2  4  6  8  10  12  14  16  18  20  22  90  158  160  162  164  166  168  170  172  174  176  178 

180]; 

cl_tab=[0.  0.  0.  0.  0.  0.  0.  0.  0.; 

.205  .205  .205  .205  .205  .205  .205  .205  .205; 


.41 

.41 

.41 

.41 

.41 

.41 

.41 

.41 

.41; 

.6 

.6 

.6 

.6  .6 

.6 

.6 

.6 

.6; 

.77 

.77 

.77 

.77 

.77 

.77 

.77 

.77 

.77; 

.82 

.82 

.82 

.82 

.82 

.82 

.82 

.82 

.82; 

.82 

.82 

.82 

.82 

.82 

.82 

.82 

.82 

.82; 

.8 

.8 

.8 

00 

00 

.8 

.8 

.8 

.8; 

.76 

.76 

.76 

.76 

.76 

.76 

.76 

.76 

.76; 

.705 

.705  .705  .705  .705  .705 

.705 

.705  .705; 

.65 

.65 

.65 

.65 

.65 

.65 

.65 

.65 

.65; 

.65 

.65 

.65 

.65 

.65 

.65 

.65 

.65 

.65; 

-.0627  -.0627  -.0627  -.0627  -.0627  -.0627  -.0627  -.0627  -.0627; 
-.98  -.98  -.98  -.914  -.934  -.926  -.875  -.838  -.822; 

-.975  -.975  -.96  -.910  -.93  -.920  -.856  -.81  -.79; 

-.969  -.969  -.962  -.906  -.926  -.914  -.838  -.782  -.758; 

-.963  -.963  -.966  -.902  -.922  -.908  -.819  -.754  -.726; 

-1.07  -1.07  -.824  -.803  -.805  -.88  -.8  -.726  -.694; 

-.718  -.718  -.532  -.528  -.66  -.83  -.79  -.698  -.662; 

-.366  -.366  -.24  -.4  -.61  -.78  -.81  -.67  -.630; 

-.245  -.245  -.3  -.33  -.55  -.74  -.75  -.666  -.622; 

-.39  -.39  -.44  -.32  -.52  -.68  -.69  -.663  -.615; 

-.4  -.4  -.42  -.44  -.47  -.58  -.47  -.486  -.428; 

-.185  -.185  -.185  -.196  -.199  -.255  -.25  -.31  -.24; 

.029  .029  .048  .048  .072  .07  .07  -.15  -.05; 

.244  .244  .282  .292  .343  .395  .35  .138  .2; 

.459  .459  .515  .536  .614  .72  .56  .39  .449; 

.673  .673  .749  .78  .84  .83  .705  .64  .7; 

.888  .888  .983  .96  .91  .882  .805  .765  .806; 

1.103  1.103  1.17  1.01  .946  .92  .842  .81  .85; 

1.25  1.25  1.13  .96  1.  .924  .845  .829  .865; 

1.1  1.1  1.03  1.07  1.053  .928  .848  .848  .88; 

.98  .98  .96  1.06  1.075  .92  .860  .867  .895; 

.982  .982  .966  1.07  1.064  .9  .880  .886  .91; 

.984  .984  .972  1.065  1.053  .9  .900  .905  .925; 

.987  .987  .979  1.06  1.042  .92  .920  .924  .94; 

.0627  .0627  .0627  .0627  .0627  .0627  .0627  .0627  .0627; 
-.66  -.66  -.66  -.66  -.66  -.66  -.66  -.66  -.66; 

-.655  -.655  -.655  -.655  -.655  -.655  -.655  -.655  -.655; 
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-.685  -.685  -.685  -.685  -.685  -.685  -.685  -.685  -.685; 
-.73  -.73  -.73  -.73  -.73  -.73  -.73  -.73  -.73; 

-.77  -.77  -.77  -.77  -.77  -.77  -.77  -.77  -.77; 


-.805  -.805  -.805  -.805  -.805  -.805  -.805  -.805  -.805; 

-.79  -.79  -.79  -.79  -.79  -.79  -.79  -.79  -.79; 

-.61  -.61  -.61  -.61  -.61  -.61  -.61  -.61  -.61; 

-.42  -.42  -.42  -.42  -.42  -.42  -.42  -.42  -.42; 

-.21  -.21  -.21  -.21  -.21  -.21  -.21  -.21  -.21; 

0.  0.  0.  0.  0.  0.  0.  0.  0.]; 

Mach_CD=[0  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0]; 

aoa_CD=[-180  -178  -176  -174  -172  -170  -168  -166  -164  -162  -160  -158  -135  -90  -60  -45  -30  -22  - 
20-18-16-14-12-10-8... 

-6  -4  -2  0  2  4  6  8  10  12  14  16  18  20  22  30  45  60  90  135  158  160  162  164  166  168  170  172 


174  176  178  180]; 


cd  tab=[. 

02  .< 

02  .02 

.02 

.02 

.02 

.02 

.02 

.02; 

.03 

.03 

.03 

.03 

.03 

.03 

.03 

.03 

.03; 

.05 

.05 

.05 

.05 

.05 

.05 

.05 

.05 

.05; 

.08 

.08 

.08 

.08 

.08 

.08 

.08 

.08 

.08; 

.11 

.11 

.11 

.11 

.11 

.11 

.11 

.11 

.11; 

.14 

.14 

.14 

.14 

.14 

.14 

.14 

.14 

.14; 

.185 

.185 

.185 

.185  .185  . 

185 

.185 

.185 

.185; 

.235 

.235 

.235 

.235  .235  .: 

235 

.235 

.235 

.235; 

.25 

.25 

.25 

.25 

.25 

.25 

.25 

.25 

.25; 

.265 

.265 

.265 

.265  .265  .: 

265 

.265 

.265 

.265; 

.295 

.295 

.295 

.295  .295  .: 

295 

.295 

.295 

.295; 

.36 

.36 

.36 

.36 

.36 

.36 

.36 

.36 

.36; 

1.1945  1.1945  1.1945  1.1945  1.1945  1.1945  1.1945  1.1945  1.1945; 
2.022  2.022  2.022  2.022  2.022  2.022  2.022  2.022  2.022; 

1.662  1.662  1.662  1.662  1.662  1.662  1.662  1.662  1.662; 

1.194  1.194  1.194  1.194  1.194  1.194  1.194  1.194  1.194; 

.6  .6  .6  .6  .6  .6  .6  .6  .6; 

.3438  .3438  .3885  .4065  .414  .458  .479  .497  .514; 

.2723  .2723  .3281  .3506  .36  .415  .441  .463  .486; 

.2007  .2007  .2678  .2948  .3267  .372  .403  .43  .457; 

.1292  .1292  .2073  .2388  .2887  .329  .3655  .397  .428; 

.0576  .0576  .147  .183  .246  .286  .3278  .363  .399; 

.0174  .0174  .0225  .12  .191  .243  .29  .33  .37; 

.008  .008  .0132  .068  .127  .177  .225  .262  .297; 

.0082  .0082  .0095  .0206  .07  .113  .16  .203  .248; 

.0079  .0079  .0085  .0097  .026  .06  .1  .149  .202; 

.0075  .0075  .008  .008  .0125  .03  .065  .115  .152; 

.0075  .0075  .008  .0075  .0085  .012  .028  .066  .117; 

.0075  .0075  .008  .0075  .008  .008  .017  .05  .09; 

.008  .008  .0082  .0075  .0075  .0105  .04  .08  .1175; 

.0085  .0085  .0085  .008  .011  .036  .09  .12  .1525; 

.009  .009  .0105  .011  .029  .081  .128  .167  .203; 

.011  .011  .014  .026  .0743  .126  .17  .21  .249; 

.017  .017  .021  .08  .1247  .162  .225  .262  .298; 

.026  .026  .0935  .153  .18  .238  .285  .3225  .363; 

.145  .145  .1635  .2121  .246  .284  .326  .357  .393; 

.2147  .2147  .2259  .2643  .2887  .329  .3655  .391  .423; 

.274  .274  .2836  .3166  .3267  .327  .403  .43  .457; 

.3333  .3333  .3414  .3688  .36  .415  .441  .463  .486; 

.2927  .2927  .3991  .421  .414  .458  .479  .497  .514; 
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.6  .6  .6  .6  .6  .6  .6  .6  .6; 

1.194  1.194  1.194  1.194  1.194  1.194  1.194  1.194  1.194; 

1.662  1.662  1.662  1.662  1.662  1.662  1.662  1.662  1.662; 

2.022  2.022  2.022  2.022  2.022  2.022  2.022  2.022  2.022; 

1.1945  1.1945  1.1945  1.1945  1.1945  1.1945  1.1945  1.1945  1.1945 


.36 

.36 

.36 

.36 

.36 

.36 

.36 

.36 

.36; 

.295 

.295 

.295 

.295  . 

295 

295 

.295 

.295 

.295; 

.265 

.265 

.265 

.265  . 

265 

265 

.265 

.265 

.265; 

.25 

.25 

.25 

.25 

.25 

.25 

.25 

.25 

.25; 

.235 

.235 

.235 

.235  . 

235 

235 

.235 

.235 

.235; 

.185 

.185 

.185 

.185  . 

185 

185 

.185 

.185 

.185; 

.14 

.14 

.14 

.14 

.14 

.14 

.14 

.14 

.14; 

.11 

.11 

.11 

.11 

.11 

.11 

.11 

.11 

.11; 

.08 

.08 

.08 

.08 

.08 

.08 

.08 

.08 

.08; 

.05 

.05 

.05 

.05 

.05 

.05 

.05 

.05 

.05; 

.03 

.03 

.03 

.03 

.03 

.03 

.03 

.03 

.03; 

.02 

.02 

.02 

.02 

.02 

.02 

.02 

.02 

.02]; 

a=alpha.*  180/pi; 
for  j  =  1  :leneth(a) 
ifaa)<-180 
aG)  =  aG)  +  360; 
elseifaG)  >180 
aG)  =  aG)  -  360; 

end 

end 

maeh  =  abs(maeh); 

CL  =  interp2(Maeh_CL,aoa_CL,el_tab,maeh,a); 
CD  =  interp2(Maeh_CD,aoa_CD,ed_tab,maeh,a); 
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APPENDIX  M.  VR12CLCD.M 


function  [CL,CD]=vrl2clcd(alpha) 

%  vrl2clcd  calculates  CL  and  CD  for  the  VR-12  airfoil 
%  given  angle  of  attack  (alpha)  in  radians 
%  [CL,CD]=vrl2clcd(alpha) 

CL=zeros(size(alpha)); 

CD=zeros(size(alpha)); 
a=alpha*  180/pi; 


chk=(a>=20  &  a<=180); 

CL=CL+chk.*(1.1733-0.018879*a+1.5752e-3*a.^2-3.1925e-5*a.^3+2.0949e-7*a.M-4.3807e- 

10*a.^5); 


chk=(a>=-180  &  a<=-50); 

CL=CL+chk.*(-4.6183-0.1923*a-3.5554e-3*a.^2-3.3273e-5*a.^3-1.4528e-7*a.M-2.3003e- 

10*a.^5); 


chk=(a>-50  &  a  <-30); 

CL=CL+chk.*(-0.22114+0.020857*a+2.8571e-4*a.^2); 


chk=(a>=-30  &  a  <=-10); 

CL=CL+chk.*(-l.ll-0.12383*a-0.01515*a.^2-6.8667e-4*a.^3-le-5*a.M); 
chk=(a<20  &  a>-10); 

CL=CL+chk.*(0.11976+0.12341*a+5.5841e-4*a.^2-2.0652e-4*a.^3); 


chk=(a>=17&a<=180); 

CD=CD+chk.*(-0.26376+0.017917*a+6.9927e-4*a.^2-9.1137e-6*a.^3+2.6277e-8*a.M); 


chk=(a>=-180  &  a<=-10); 

CD=CD+chk.*(-0.17486-0.034463*a-1.0233e-4*a.^2-2.8958e-6*a.^3-4.6577e-8*a.M-1.5557e- 

10*a.^5); 


chk=(a>-10  &  a<=0); 

CD=CD+chk.*(9.8678e-3+3.4934e-3*a+1.4844e-3*a.^2-1.3564e-4*a.^3-1.0936e-5*a.M); 


chk=(a>0  &  a<=15); 

CD=CD+chk.*(9.8e-3+7.0457e-4*a+5.6104e-5*a.^2-4.1151e-5*a.^3+3.8695e-6*a.M); 
chk=(a>15  &  a<17); 

CD=CD+chk.*(-1.33+1.325e-l*a-2.5e-3*a.^2); 
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